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ABSTRACT 

During  the  lifetime  of  this  project,  the  research  team  significantly  advanced  the  state  of  understanding  of  modu¬ 
lated  polarimeter  systems  for  active  and  passive,  imaging-  and  non-imaging  polarimetry.  Modulated  polarimeters 
infer  the  polarization  state  of  light  by  creating  a  set  of  polarization-dependent  carriers  that  are  modulated  by 
the  intensity  signal.  These  carriers  can  be  created  in  any  independent  domain,  such  as  time,  space,  wavenumber, 
angle  of  incidence,  and  they  can  be  created  in  combinations  of  domains  simultaneously.  The  work  supported  on 
this  project  has  solidified  the  theory  behind  such  instruments,  allowing  for  new  design  philosophies  that  improve 
state-of-the-art  instruments.  When  the  project  began  five  years  ago,  only  a  cursory  understanding  of  modulated 
instruments  existed,  and  the  data  reduction  matrix  was  the  primary  means  of  processing  polarization  data.  In 
this  project,  the  DRM  has  been  expanded  to  include  a  full  functional  formalism,  allowing  for  a  range  of  new 
tools  in  polarimeter  design  to  be  brought  to  bear.  During  the  course  of  this  project,  five  PhD  students,  seven  MS 
students,  and  two  undergraduates  were  trained.  Four  PhD  dissertations  were  primarily  supported  by  funding 
from  this  project. 


1.  EXECUTIVE  SUMMARY 
1.1  Principal  Accomplishments 

This  project  was  carried  out  in  the  Advanced  Sensing  Laboratory  of  the  College  of  Optical  Sciences  at  the 
University  of  Arizona.  The  original  project  dates  were  April  15,  2010  -  December  15,  2015.  The  was  initiated 
based  on  work  that  the  PI  had  done  under  previous  USAF  funding  with  microgrid  imaging  polarimeters 
While  the  main  motivation  at  the  beginning  of  the  project  was  to  improve  the  performance  of  microgrid  instru¬ 
ments,  the  project  expanded  signficantly  to  include  all  forms  of  modulated  polarimeters  such  as  time  modulated 
instruments,®  prismatic  polarimeters,®  spectrally  channeled  polarimeters,^  and  even  instruments  modulated  in 
multiple  domains 00  It  is  important  to  note  that  much  of  the  progress  was  achieved  in  collaboration  with  Pro¬ 
fessor  Russell  Chipman  of  the  University  of  Arizona.  Prof.  Chipman  did  not  receive  any  direct  funding,  from 
this  project,  but  he  co-supervised  one  of  the  students  and  was  an  active  collaborator  on  several  of  the  projects. 

The  progress  supported  on  this  grant  can  be  broken  down  into  the  following  areas: 

•  Fundamental  understanding  of  modulated  instruments:  Dr.  Charlie  LaCasse  worked  on  a  basic  theory  of 
modoluated  instruments  for  his  PhD  dissertation.  That  work  extended  our  earlier  work  with  microgrids  and 
coupled  it  with  Prof.  Chipman’s  work  on  temporally  modulated  instruments0E2l  Dr.  LaCasse  extended 
this  work  to  cover  random  signal  processing  and  the  associated  transfer  function  formalisms.11 

•  Scene-Based  Nonuniformity  Correction  for  Focal  Plane  Arrays:  Dr.  Wiley  Black  worked  on  methods  to 
correct  residual  calibration  errors  in  the  FPA  due  to  nonuniformity  in  the  array  response.  Much  of  his  work 
turned  out  to  be  generally  applicable  to  all  imagers Ji2  13  but  he  also  considered  its  specific  application  to 

infrared  microgrid  instrument^® 
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•  Coherence  Manipulation  for  Sensing-.  Dr.  Oscar  Rodriguez-Herrera  developed  strategies  to  be  able  to 
manipulate  the  coherence  of  light  in  order  to  indirectly  sense  BRDF  from  monostatic  measurements. 10 
This  was  follow-on  work  from  preliminary  results  obtained  under  earlier  AFOSR  funding,16  and  continues 
in  current  work. 

•  Generalized  Channeled  Polarimetry :  Dr.  Andrey  Alenin  developed  a  general  theoretical  framework  that 
allows  channeled  polarimeters  to  be  designed  from  a  priori  information.1'  Previous  design  methods  were 
ad  hoc,  and  generally  resulted  in  needlessly  sub-optimal  systems.  This  work  was  a  direct  outgrowth  of  Dr. 
LaCasse’s  research  (mentioned  above). 

•  Partial  Mueller  Matrix  Polarimeters  (pMMPs):  Dr.  Alenin  also  worked  on  the  development  of  pMMPs.18 
These  devices  are  active  instruments  that  intentionally  omit  certain  Mueller  elements  in  order  to  improve 
bandwidth  and  were  pioneered  by  the  PI  and  Brian  Hoover.19  20  He  also  developed  methods  to  apply  his 
generalized  channeled  polarimeter  theory  to  pMMPsP1  Results  from  this  work  were  implemented  in  an 
instrument  built  at  AFRL/RYJT  by  a  student  from  the  Pi’s  group.22 

•  Hybrid  Modulated  Polarimeters:  As  part  of  his  research,  Dr.  LaCasse  explored  polarimeters  modulated  in 
both  space  and  time. '  Dr.  Israel  Vaughn  extended  this  work  to  include  the  design  and  optimization  of  active 
spatiotemporally  modulated  systems, and  he  built  a  portable  polarimeter  based  on  these  principles.^ 
These  techniques  have  also  been  extended  by  us  to  hybrid  wavelength  and  time  modulated  instruments  in 
collaboration  with  Frans  Snik  and  colleagues  in  the  Netherlands  for  astronomical  applications. El 

•  Other  Collaborative  Works:  The  PI  and  Drs.  Vaughn  and  Rodriguez-Herrera  had  an  ongoing  collaboration 
with  Prof.  Toshitaka  Wakayama  from  Saitama  Medical  University  in  Japan.  The  collaboration  began  in 
2012  when  Dr.  Wakayama  visited  the  laboratory  for  three  months.  Our  group  worked  with  Dr.  Wakayama 
to  develop  axial  waveplates  in  the  terahertz20  and  to  apply  these  instruments  to  angle-modulated  sensing.26 
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2.  INTRODUCTION 

Optical  polarization  sensing  is  an  example  of  an  indirect  imaging  or  sensing  method,^  since  optical  detectors 
are  generally  insensitive  to  the  polarization  state  of  the  incident  radiation.  Because  of  that,  the  nature  of  the 
performance  of  polarimeters  -  including  SNR,  accuracy,  calibration  error,  etc.  -  is  intrinsically  related  to  the 
processing  methods  used  to  compute  the  desired  information  from  the  actual  measured  intensity  signals. 

A  widely  used  general  formalism  for  describing  the  polarization  state  of  optical  radiation  in  the  partially 
coherent  case  is  the  Stokes-Mueller  calculus.^  in  this  formalism,  the  polarization  state  of  light  is  described  by 
a  set  of  four  Stokes  parameters  that  are  related  to  the  second  statistical  moments  of  the  field.  These  moments 
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can  be  related  to  any  radiometric  quantity,  but  for  imaging  polarimeters,  the  irradiance  is  a  reasonable  quantity 
and  the  Stokes  parameters  can  be  defined  as 


S  = 


lo  + 19  0 
lo  ~  h  0 

I‘ 15  —  -^135 

I L  -  I R 


(2.1) 


where  Ix,  Iy ,  J45,  and  /135  are  the  observed  irradiances  through  linear  polarizers  oriented  at  0°,  45°,  90°,  and 
135°,  respectively,  and  I ^  and  Ir  are  the  irradiances  observed  through  left-  and  right-circular  polarizers.  Note 
that  the  angles  are  defined  with  respect  to  an  arbitrary  coordinate  system.  When  light  interacts  with  an  object 
or  material,  the  polarization  properties  are  altered.  When  the  interactions  are  linear,  the  changes  in  the  Stokes 
parameters  can  be  described  by  a  4  x  4  real  matrix  called  the  Mueller  matrix 
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The  notation  y  0  J  indicates  that  the  polarimetric  quantities  are  functions  of  some  set  of  independent  parameters 
that  may  include  spatial  coordinates,  time,  wavenumber,  angle  of  incidence,  etc. 


2.1  Wavefront  Division  Polarimeters 

A  Stokes  polarimeter  is  used  to  estimate  a  set  of  Stokes  parameters.  There  are  two  broad  classes  of  polarimeters: 
wavefront  division  polarimeters  and  modulated  polarimeters.  Wavefront  division  polarimeters  operate  by  making 
multiple  copies  of  the  incident  wave  field  and  relaying  them  to  separate  detectors  through  different  polarization 
optics.  The  most  common  examples  of  wavefront  division  polarimeters  are  Division  of  Amplitude  polarimeters, 
which  use  beam  splitters  to  create  the  copies  of  the  field,  and  Division  of  Aperture  Polarimeters,  which  use  an 
array  of  optics  to  split  the  wavefront  in  an  aperture  planePS  Subsequet  to  this  splitting,  both  strategies  employ 
independent  polarization  optics  in  each  of  the  channels  in  order  to  affect  the  measurement.  Each  channel  is 
described  by  a  polarization  analyzer  vector  A.,,  and  the  power  measured  in  the  channel  is 


Xi  =  A  ■  S. 


(2.3) 


The  set  of  measurements  can  be  compiled  into  a  linear  system  of  equations 


=  W  S, 


(2.4) 


and  the  Stokes  parameters  are  estimated  as 

S  =  W+  •  X.  (2.5) 

The  matrix  W+  is  a  suitable  inverse  or  pseudoinverse,  and  is  often  referred  to  as  the  data  reduction  matrix.'28 


2.2  Modulated  Polarimeters 

A  wavefront  division  polarimeter  splits  the  light  into  channels  that  are  measured  by  different  detectors  (or  subsets 
of  detector  elements  in  some  DoAP  systemJ^).  I11  contrast,  a  modulated  polarimeter  uses  a  single  detector  or 
array  of  detectors  to  make  all  measurements,  but  uses  polarization-dependent  carriers  to  create  channels  that 
carry  the  polarization  information.  Most  systems  have  employed  purely  periodic  carriers,  which  results  in  an 
orthogonal  frequency  division  multiplexed  system,  though  non-periodic  carriers  can  be  employed  as  well.8 

The  detector  output  of  a  modulated  polarimeter  (which  is  proportional  to  irradiance  or  power  received)  is 
X{t)  =  A(t)  ■  S (t)  =  Ao(t)so(t)  +  A1(t)s1(t)  +  A2(t)s2(t)  +  A3(t)s3(t).  (2.6) 


Iii  Eq.  2.6  we  have  only  indicated  temporal  modulation,  or  ^  9  ^  =  ( t ).  However,  the  theory  is  immediately 
applicable  to  any  combination  of  modulation  domains.  When  the  analyzer  vector  is  periodic,  then  the  Fourier 
transform  of  Eq.  2.6  produces  a  set  of  channels  that  are  orthogonal  in  the  Fourier  domain.  The  Stokes  parameters 
can  then  be  estimated  through  a  set  of  demodulation  and  filtering  operations. 


2.3  Mueller  Polarimeters 

It  should  be  noted  that  the  Stokes  polarimeter  can  be  generalized  into  an  active  Mueller  polarimeter.31;  This 
instrument  includes  both  a  polarization  state  analyzer  (PSA),  which  is  a  Stokes  polarimeter,  and  an  polarization 
state  generator  (PSG).  The  PSG  is  a  PSA  operated  in  reverse  that  generates  a  controlled  illumination  state. 
The  intensity  measured  for  a  particular  PSG/PSA  combination  is 

Xi  =  =  Dj  ■  M.  (2.7) 

The  vector  Dj  is  a  16  x  1  vector  obtained  by  rearranging  the  elements  of  the  dyad  product  GA1,  and  M  is 
obtained  by  reordering  the  elements  of  M.  All  of  the  analyses  that  follow  are  equally  applicable  to  both  Stokes 
and  Mueller  polarimeters. 


3.  MODULATED  POLARIMETER  DATA  PROCESSING 
3.1  Data  Reduction  Matrix  method  for  DoT  Polarimeters 

The  conventional  theory  of  operation  for  DoT  polarimeters  (and  for  other  classes  of  polarimeter  as  well)  is  well 
known—  and  is  briefly  reproduced  here.  The  polarimeter  operates  by  using  a  collection  of  retarders  and  analyzers 
that  can  be  changed  from  measurement  to  measurement.  The  most  common  example  is  the  rotation  of  a  wave 
plate  in  front  of  a  linear  analyzer,^2  but  other  optical  configurations  can  also  be  used.33lMl  The  incident  Stokes 
vector  is  altered  by  the  Mueller  matrix  of  the  polarization  optics,  and  the  modified  irradiance  is  measured  by  a 
polarization-blind  detector.  This  irradiance  is  given  mathematically  by  taking  the  inner  product  of  the  first  row 
of  the  system  Mueller  matrix  with  the  incident  Stokes  vector  as 

m  =  sA(t)Tsin,  (3.i) 

where  SA  (t)  is  the  time- varying  anaylzer  Stokes  vector  given  by  the  first  row  of  the  system  Mueller  matrix.  For 
the  purposes  of  the  present  discussion,  we  will  assume  that  the  irradiance  is  sampled  ideally.  The  effects  of  an 
integration  time  that  is  comparable  to  the  time  variation  of  S A(t)  can  be  included  using  standard  methods.  The 
nth  sample  of  the  irradiance  is  then 

I[n]  =  SA[n]TS  =  sao[™]so  +  SAi[n]si  +  sA2[n]s2  +  SA3[«]s3-  (3.2) 


In  general  practice,  it  is  assumed  that  the  polarimeter  modulates  the  irradiance  at  rates  that  are  much  faster 
than  the  temporal  variation  of  the  incident  Stokes  vectorJH  35  and  that  the  Stokes  parameters  in  Eq. 


6.1 


constants.  When  this  is  the  case,  we  can  take  a  series  of  measurements  and  form  a  system  of  linear  equations 

(3.3) 


I  = 


'  m 

'  Sa[1]t  ' 

.  m . 

.  sA[N}T  _ 

Sm  =  WS; 


The  incident  Stokes  vector  is  determined  by  matrix  inversion  (or  pseudoinversion)  as 


sm  =  ^-T=(^TwfVi. 


(3.4) 


where  the  hat  in  Eq.  3.4  indicates  that  the  result  is  an  estimate  of  the  true  value.  The  matrix  W  is  generally 
referred  to  as  the  data  reduction  matrix  (DRM).28 


While  the  above  discussion  approaches  the  problem  from  a  linear  algebra  perspective,  a  frequency  domain 
approach  is  often  more  informative.  Equation  |3.3|is  easiest  to  understand  with  the  simple  example  of  the  rotating 
analyzer  linear  polarimeter.  If  we  rotate  an  ideal  linear  analyzer  at  a  constant  angular  frequency  ujq  (in  units  of 
degrees  per  sample)  we  have 


and  Eq.  |6.1|  becomes 


sA[n\  =  \[  1 


=  2  ts° 


cos(2oj0n)  sin(2w0^)  0 


cos(2wo«)si  +  sin(2wo^)s2]  ■ 


(3-5) 


(3-6) 


Equation  |3.6| tells  us  that  the  irradiance  information  is  carried  in  the  DC  term,  the  Si  information  is  carried 
in  the  in-phase  2wo  term,  and  the  S2  information  is  carried  in  the  quadrature  2wo  term.  Since  this  polarimeter 
does  not  have  a  retardance  element,  it  cannot  measure  S3.  The  polarization  information  could  be  demodulated 
using  frequency-domain  methods,  which  is  done  for  some  systems  using  photo-elastic  modulators  operating  at 
high  rates  (~10s  of  kHz)P^  However,  the  linear  algebra  formulation  of  Eq.  3.4  is  often  preferred  since  it  is  easier 


to  incorporate  calibration,  error  analysis,  and  other  pratical  considerations.  In  principle,  only  four  measurements 
are  required  to  uniquely  invert  Eq.  |3.4|  However,  it  is  common  to  employ  more  than  4  measurements  in  order 
to  help  improve  the  accuracy  of  the  polarimeter  and  reduce  the  effects  of  noise  and  systematic  errors.4  US!1  ]HD 

3.2  Using  a  Linear  Systems  Formalism  for  Processing 

With  only  a  few  exceptions,4  39  most  papers  in  the  polarimetry  literature  treat  S in  as  approximately  constant  in 
computing  the  DRM.  Diner,  et  al.,  initially  make  this  assumption  in  their  dual-PEM  polarimeter.  Their  device 
makes  40  time-sequential  measurements  to  reconstruct  So,  Si  and  S2,  resulting  in  a  highly  underdetermined 
system  of  equations  with  17  unconstrained  degrees  of  freedom.  After  initially  assuming  that  the  measured 
Stokes  vector  was  constant,  they  then  used  one  of  their  degrees  of  freedom  to  innoculate  the  DRM  to  linear 
gradients  in  the  underlying  Stokes  vector. 

An  examination  of  Eq.  |6.1|or  Eq.  |3.6|reveals  that  the  constant  Stokes  vector  assumption  is  overly  limiting,  and 
even  the  linearly-varying-in-time  assumption  of  Diner,  et  al.,  is  still  unnecessarily  limiting.  It  is  apparent  that 
basic  communications  systems  theory  can  be  employed  allowing  the  Stokes  parameters  to  be  functions  of  time  as 
well,  subject  to  a  bandwidth  criterion  imposed  by  the  reconstruction  process.  This  result  is  analagous  to  earlier 
developments  for  spatially  modulated  imaging  polarimeter£E2l  and  spectrally  modulated  spectropolarimetersP^ 


To  introduce  this  concept  for  DoT  polarimeters,  we  take  the  Fourier  transform  of  the  signal  I(t)  in  Eq.  3.1 
but  now  we  allow  the  Stokes  parameters  to  become  functions  of  time  as 

/(w)  =  §ao{u)  *  s0(w)  +  SAl(w)  *  Si(w)  +  SA2(u)  *  S2(w)  +  §A 3(w)  *  s3(w),  (3.7) 

where  the  tilde  indicates  the  Fourier  transform  and  *  is  the  convolution  operator. 

Equation  |3.7|  presents  a  conventional  deconvolution  problem  that  can  be  inverted  by  careful  design  of  the 


analyzer  Stokes  vector  SA(t).  The  specific  modulation  strategy  of  Eq.  3.6  creates  distinct  side  bands  in  frequency 


space  that  carry  the  polarization  information.  Using  the  expressions  of  Eq.  |3.6|in  Eq.  |3.7| produces 
I(uj)  =  —  ^so(w)  +  —  (si(w  —  2wq)  +  Si(w  +  2ojq))  +  —  (s2(w  —  2wq)  —  S2  (u>  +  2wq)) 


(3-8) 


We  see  from  Eq.  3.8  that  forcing  the  Stokes  vector  to  be  constant  is  a  severe,  and  unnecessary,  restriction.  If  we 
assume  that  so  is  band  limited  to  frequencies  below  Wo  and  si  and  s2  are  band  limited  to  frequencies  below  Wi, 
we  can  accomplish  error-free  reconstruction  (in  the  noise-free  case)  provided  that 


W0  +  Wi  <  wa/2, 


(3-9) 


where  ws/27t  is  the  temporal  sampling  frequency  of  the  polarimeter.  Equation  3.9  is  the  Nyquist  sampling 


criterion  for  a  rotating  analyzer  polarimeter.  It  is  worth  noting  that  the  band  limit  requirement  is  placed  on 
Wo  and  Wi  together,  so  it  is  possible  to  apportion  a  given  sampling  bandwidth  among  the  various  channels 
in  question.  Other  types  of  DoT  polarimeters  with  different  specific  modulations  will  have  band  limit  criteria 
modified  somewhat  from  Eq.  |3.9|  but  similar  critera  can  be  derived  for  all  modulated  polarimeters. 


3.3  Unifying  the  DRM  and  linear  systems  methods 

In  this  section  we  generalize  the  linear  systems  formalism  to  encompass  the  DRM  method.  Doing  this  shows 
an  immediate  drawback  of  the  DRM  method  that  leads  to  unncessary  polarimetric  aliasing  that  can  easily  be 
remedied  by  choosing  a  different  weighting  scheme  in  formulating  the  DRM. 

Consider  a  Stokes  vector  signal  to  be  measured  that  is  a  function  of  space,  time,  and  wavelength  S in(x,  y, t,  A). 
The  system  has  an  impulse  response  function  that  is  described  by  h(x,  y,  t,  A),  a  detector  with  a  sampling  function 
d(x,y,t,  A),  and  a  polarimetric  modultation  described  by  SA(x,  y,  t,  A).  We  also  consider  a  uniform  sampling 
process  in  space,  time,  and  wavelength  with  sampling  intervals  X,  Y .  T,  and  43  Nonuniform  sampling  is 
possible  as  well,  but  complicates  the  analysis  unnecessarily.  With  these  definitions,  the  modulated  irradiance  is 
given  as 


I(x,y,t,X)  =  [  [SA(x,  y,  t,  A)T  (h(x,y,t,  A)  *  Sin(x,  y,  t,  A))]  *  d(x,y,t,  A)]  comb 


f  x_  y_  t_  A  A 

\X’  Y’  T:  L  J 


(3.10) 


In  Eq.  |3.10|the  operator  *  defines  convolution;  however,  the  various  functions  are  not  necessarily  linear,  shift- 
invariant  (LSI)  functions  in  general.  When  the  system  is  not  LSI,  the  convolution  integrals  have  to  take  the 
more  general  form 

/OO 

f(a)h(x,a)da  (3-11) 

-CO 

rather  than  the  more  familiar  form  for  LSI  systems 


92{x) 


f{a)h{x 


a)  da, 


(3.12) 


where  in  both  cases  the  input  is  f{x)  and  the  impulse  response  is  h(x).  Furthermore,  the  point  spread  function 
h(x,y,t,  A)  is  assumed  to  be  scalar.  In  reality,  the  system  is  really  described  by  a  polarimetric  impulse  response 
matrix  that  describes  how  the  optical  system  alters  the  polarization  state  between  object  and  image  plane  before 
it  is  ever  sampled.^  For  convenience  we  will  define  our  units  so  that  X  =  Y  =  1,  T  =  1,  and  L  =  1,  so  we  will 
be  working  in  normalized  frequency  space  (cycles  per  sample). 

Referring  back  to  Eq.  |3.4|  the  DRM  formalism  computes  the  pseudoinverse  as 

yg  1  =  (^T^)  (3.13) 


and  the  irradiance  I{x,y,t,  A)  is  rearranged  into  a  column  vector  that  W  1  operates  on.  The  first  part  of  the 
pseudoinverse  in  Eq.  |3.13|  is 

1=  (^T^)_1=  ^SA(n)SA(n)Tj  .  (3.14) 

In  this  notation,  the  polarimeter  uses  N  measurements  to  form  the  measurement  matrix.  For  example,  a  rotating 
retarder  polarimeter  might  make  TV  =  16  measurements  as  the  retarder  is  rotated  from  0°  to  360°,  or  a  DoFP 
polarimeter  might  be  decomposed  into  4-element  (2  x  2)  superpixels.  In  the  linear  systems  formalism,  we  propose 
a  similar  formation  of  Z 

^=1  ^3  sA{x,y,t,X)SA(x,y,t,  A)T  j  .  (3.15) 

\x,y,t,  A  J 

The  WTW  part  of  the  inverse  is  the  integral  of  the  product  of  the  modulation  functions  contained  in  SA,  so  the 
inversion  of  this  quantity  describes  how  to  separate  the  polarization  signals  if  their  modulation  is  not  completely 
orthogonal.  For  modulation  schemes  that  are  orthogonal  over  the  integral  the  quantity  WTW  and  its  inversion 
will  simply  be  a  diagonal  matrix.  We  refer  to  Z  as  The  modulator  inner  product  inversion  matrix. 


The  system  could  equivalently  be  uniformly  sampled  in  wavenumber  1/A,  which  is  also  common. 


If  the  measurement  has  been  taken  such  that  there  is  an  integer  number  of  periods  of  modulation  in  all  of 
the  functions  contained  in  a,  this  quantity  will  be  a  constant  with  respect  to  the  initial  phase  of  the  modulation. 
However,  in  cases  where  the  polarimeter  varies  in  time  or  in  space,  as  would  be  the  case  when  there  is  drift  in  the 
absolute  angular  position  of  a  rotating  retarder  or  in  the  absolute  phase  of  oscillation  of  the  dual-PEM  system,4 
then  it  is  possible  that  Z  might  need  to  be  computed  locally.^ 


Next  we  turn  to  the  W1  term  in  Eq. 


3.14 


To  understand  this  element  it  is  useful  to  examine  how  this  term 


operates  on  the  modulated  irradiance  in  the  standard  DRM  formalism.  We  have 


WTI  =  WTWSn. 


(3.16) 


Examining  this  equation  closely  in  the  context  of  the  linear  systems  formalism  reveals  that  the  matrix  WT  plays 
two  roles.  The  first  role  is  the  homodyne  in  the  demodulation  process.  The  matrix  W  includes  the  modulation 
strategy  of  SA(x,y,t,  X),  multiplying  the  input  signal  by  W  then  WT  is  equivalent  to  mixing  with  a  carrier 
frequency  once  to  move  the  base  band  signal  up  to  the  side  bands  and  a  second  time  to  create  a  copy  at  base 
band  along  with  spurious  copies  at  higher  frequencies.  In  communications  theory,  the  next  step  is  to  low  pass 
filter.  This  low  pass  filter  is  implicitly  included  in  the  matrix  multiplication  WTW.  However,  in  the  matrix 
multiplication  the  low  pass  filter  has  a  rectangular  footprint  (in  time,  space,  wavelength).  This  rectangular 
footprint  results  in  a  sing-function  low  pass  filter  that  results  in  leakage  of  spurious  high  frequency  signals  into 
the  reconstruction  as  we  will  see  below. 


In  the  linear  systems  formalism  we  want  to  leave  ourselves  the  freedom  to  use  an  arbitrary  low-pass  filter. 
Specifically,  we  want  to  use  band-limited  low  pass  filters  that  are  designed  to  only  pass  the  base  band  signal, 
thus  eliminating  polarimetric  aliasing.  For  the  linear  systems  formalism,  we  accomplish  this  by  separating  the 
homodyne  from  the  low  pass  filter  as 

WTI  -»  m(x,  y,  t,  A)  *  SA(x,y,t,X)I(x,y,t,X)  =  w(x,y,t,  A)  *  SA(x,  y,  t,  X)SA(x,  y,  t,  A)TSm,  (3.17) 

where  w(x,  y,  t,  A)  is  the  windowing  function  in  the  space-time- wavelength  domain  that  corresponds  to  the  desired 
low  pass  filter  in  the  corresponding  frequency  domain. 

Putting  all  of  the  terms  together,  we  can  form  the  linear  systems  version  of  the  DRM  method.  The  estimated 
Stokes  vector  distribution  is 


S  =  W  {I(x,y,t,  X)}  =  w(x,  y,  t,  A)  *  ZSA(x,y,t,  X)I{x,y,t,  X). 


(3.18) 


3.4  Discussion 

In  order  to  understand  the  polarimetric  demodulation  operation  from  the  linear  systems  perspective,  we  consider 
the  two  examples  of  a  rotating  retarder  polarimeter  and  a  dual-PEM  polarimeter. 

3.4.1  Rotating  Retarder  Polarimeter 

Consider  once  again  the  rotating  retarder  (RR)  polarimeter.  We  assume  that  the  polarimeter  is  an  ideal  linear 
retarder  of  retardance  8  followed  by  an  analyzer.  Since  our  reconstruction  was  designed  to  be  insensitive  to 
absolute  phase,  we  take  the  phase  to  be  0  with  respect  to  t  =  0.  In  this  case  the  modulation  term  for  the  kth 
temporal  sample  is 


S  A[n,m,k,l\  =  S+fc]  =  \ 
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(3.19) 


where  there  are  K  samples  per  360°  rotation  of  the  retarder.  The  modulator  inner  product  inversion  matrix  is 
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We  begin  by  assuming  that  the  excitation  has  the  form  sine  ()2  (t)  in  only  one  of  the  four  Stokes  parameters. 
Even  though  this  excitation  is  non-physical  because  it  only  includes  one  Stokes  parameter,  we  can  decompose 
an  arbitrary  physical  excitation  into  these  components  in  order  to  understand  where  the  information  from  each 


Stokes  parameter  lies  in  the  side  bands.  Figure 
parameter. 


3.1 


shows  the  measured  irradiance  for  the  sine  ()  signal  in  each 


(a) 


(b) 


(c) 


(d) 


Figure  3.1:  Shows  the  components  of  the  irradiance  due  to  S'o(a),  S'i(b),  ^2(0),  and  Sh(d) 


Figure  |3.2| shows  the  measured  irradiance  and  its  Fourier  transform  for  the  fully  polarized  signal 


S  =  sine  (t)2  [  y/3  1  1  1  ] 1  • 


(3.21) 


(a)  (b) 

Figure  3.2:  (a)  shows  the  measured  irradiance  for  S  =  [\/3111]  in  the  temporal  domain. and  (b)  shows 
the  Fourier  Transform  of  the  measured  irradiance.  Need  to  explain  colors,  and  real/imag 


As  expected,  the  irradiance  signal  shown  in  Fig.  |3.2fr  is  now  all  non-negative.  Figure  |3.2|j  shows  the  Fourier 
transform  of  the  measured  irradiance  decomposed  into  the  components  that  are  due  to  each  of  the  signals  from 

Fig-  mu 

The  first  step  in  the  reconstruction  process  it  to  multiply  the  measured  irradiance  with  the  functions  given 
in  Eq.  |3.19|in  order  to  move  the  information  for  the  desired  Stokes  parameters  down  to  base  band. 


(c) 


(d) 


Figure  3.3:  Pre  modulation  of  the  measured  irradiance  in  the  Fourier  Domain  by  ao(a),  ai(b),  02(c), 
03(d) 


Figure|373|shows  the  Fourier  transform  of  the  modulated  irradiance  using  each  of  the  four  modulation  functions 
in  Eq.  |3.19[  As  we  can  see,  the  first  and  second  modulation  functions  move  the  so  and  si  information  to  base 
band.  The  third  row  moves  the  S2  information  to  base  band,  while  the  fourth  row  moves  the  S3  information 
to  base  band.  In  order  to  separate  So  and  Si,  it  is  still  necessary  to  multiply  by  the  modulator  inner  product 
inversion  matrix  Z,  and  the  result  of  that  multiplication  is  shown  in  Fig.  3.4 


(c) 


(d) 


Figure  3.4:  Pre  filtered  Sr  for  S'ro(a),  5ri(b),  Sr2(c),  5r3(d) 


At  this  point  we  can  examine  Fig.|3.4|and  see  that  the  band  limited  information  has  been  moved  to  base  band 
independently  for  each  of  the  four  Stokes  paramters.  Now  the  only  step  left  is  to  select  a  suitable  low-pass  filter 
to  extract  the  data.  This  is  the  step  where  virtually  all  division  of  time  polarimeter  reconstruction  strategies 
make  the  crucial  mistake  of  using  a  time-limited  reconstruction  window  w(x,y,t,  A)  rather  than  a  band-limited 
one.  Figure  3.5  shows  the  effect  of  this  when  using  a  standard,  16-element  rectangular  window  for  w(x,y,t,  A). 
Each  of  the  four  panels  shows  the  Fourier  transform  of  one  of  the  Stokes  paramters  decomposed  into  the  portions 


that  arise  from  each  of  the  inputs  in  Fig.  3.1 
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(c)  a  (d)  a 

Figure  3.5:  Result  of  multiplication  of  a  filter  from  a  16  element  rect  window  in  the  Fouerier  Domain. 
Sro(a),  Srl( b),  Sr2( c),  5.3(d) 


We  can  see  that  the  proper  signal  is  reconstructed  at  base  band  in  each  case,  but  there  is  also  a  form  of 
polarimetric  aliasing  error.  This  error  manifests  both  as  self-error  (e.g.,  baseband  so  information  aliasing  up 
to  higher  frequencies  in  so  through  the  demodulation  process)  and  cross  error  (e.g.,  base  band  So  information 
showing  up  as  high  frequency  error  in  the  sj  signal).  Figure [T6] shows  the  reconstructed  Stokes  parameters  in 
the  time  domain,  where  the  high  frequency  reconstruction  error  is  evident. 


(a)  (b) 


(c)  (d) 

Figure  3.6:  Result  of  multiplication  of  a  filter  from  a  16  element  rect  window  in  the  temporal  domain. 

Sr o(a),  Sr i(b),  Sr2{ c),  s,3(d) 


3.4.2  Performance  Issues 


The  theory  presented  in  section  [373]  along  with  the  example  of  the  RR  polarimeter  in  section  |3.4.1 1  demonstrate 
some  issues  associate  with  traditional  DRM  processing  of  modulated  polarimeter  data.  The  conventional  DRM 
formulation  is  created  assuming  that  the  input  Stokes  vector  is  constant  in  the  modulated  dimension.  For 
that  reason  a  standard  matrix  formalism  is  obvious  that  equally  weights  all  of  the  observations  in  an  IV-element 
window  in  performing  the  inversion.  Since  the  excitation  in  this  case  is  assumed  DC,  the  analysis  only  guarantees 
proper  reconstruction  for  DC  signals.  In  this  paper,  we  have  demonstrated  that  there  is  a  fundamental  band 
limit  criterion  that  must  be  satisfied  for  a  polarimeter  that  is  related  to  both  the  physical  sampling  frequency 
and  the  modulation  frequency  that  introduces  the  polarization  side  bands.  So  long  as  Eq.  |3.9|is  satisfied,  error- 
free  polarimetric  reconstruction  is  possible.  This  result  relaxes  the  requirement  for  extremely  high  modulation 
frequencies  that  has  been  used  recently  to  improve  error  in  DoT  polarimeters  M 


The  analysis  of  Section |h3]used  an  analysis  of  the  DRM  formalism  to  illuminate  the  linear  systems  viewpoint. 
However,  it  is  equally  possible  to  to  use  the  knowledge  that  band  limited  reconstructions  exists  to  create  a  DRM 
picture  that  can  be  used  in  the  traditional  linear  algebra  sense.  The  two  viewpoints  must  be  identical,  but  some 
readers  might  find  one  easier  to  implement  than  the  other  for  their  purposes. 


In  order  to  construct  band-limited  DRMs  that  can  reconstruct  signals  that  obey  Eq.  3.9  it  is  necessary  to 
augment  the  rows  of  the  instrument  matrix  with  vectors  that  represent  the  Fourier  components  in  the  base  band 
used  to  reconstruct  the  Stokes  parameters. 


The  theory  presented  here  was  focused  on  DoT  polarimeters.  However,  this  scheme  applies  equally  to  all 
forms  of  modulated  polarimeters.  Similar  linear  systems  analyses  have  been  performed  for  channeled  spectropo- 
larimetersPIM  but  this  paper  represents  the  first  attempt  to  unify  the  two  theories.  Furthermore,  in  this  and 


other  analyses,  the  modulation  is  always  assumed  to  be  in  only  one  domain,  i.e.  temporal  modulation  in  a  DoT 
polarimeter,  spatial  modulation  in  a  DoFP  polarimeter,  or  spectral  modulation  in  a  channeled  spectropolarime- 
ter.  Equation  3.10  anticipates  the  case  where  the  polarimetric  modulation  represented  by  S^a:,  y,  t,  A)  actually 
modulates  across  dimensions  in  a  way  that  allows  more  accurate  reconstruction  and/or  greater  bandwidth  in 
particular  applications. 


4.  SPECTRAL  DENSITY  RESPONSE  FUNCTIONS  FOR  MODULATED 

POLARIMETERS 

4.1  Background 

Conventional  imaging  systems  are  often  described  by  their  point  spread  function  (PSF)  or  optical  transfer  function 
(OTF)  in  space  and  their  temporal  impulse  response.  These  functions  can  be  used  to  objectively  compare  systems 
based  on  criteria  that  are  determined  by  the  specific  imaging  problem.  When  comparing  modulated  polarimeters, 
transfer  functions  in  their  conventional  form  cannot  be  employed.  This  is  because  modulated  polarimeters 
multiplex  signals  in  the  frequency  domain;  a  single  frequency  sinusoidal  object  will  result  in  an  image  with 
multiple  frequencies.  This  ruins  the  linear  shift  invariant  assumption  required  for  the  transfer  function  formalism 
to  be  valid,  since  a  measurement  at  a  given  frequency  could  have  come  from  a  multiplicity  of  input  signals.  A 
response  function  for  modulated  polarimeters  is  developed  here  to  perform  an  analogous  objective  comparison 
for  polarimetric  imaging  tasks  that  the  OTF/PSF  formalism  allows  for  conventional  imaging  systems. 

The  band-limited  filter  algorithm  presented  in  section  [3]  is  to  establish  a  method  to  objectively  compare 
modulation  schemes  (spatial,  temporal,  or  some  combination)  to  determine  the  method  that  performs  “best” 
on  a  polarized  object  that  has  finite  bandwidth  both  in  the  temporal  and  spatial  domain.  This  tool  is  termed 
the  spectral  density  response  function.  The  tool  is  used  in  Section  |4.3|to  determine  an  object-specific,  optimal 
Wiener  filter  as  an  example  of  the  application  of  the  theory. 

Power  spectral  densities  (PSDs)  are  used  to  describe  objects  and  systems  when  deterministic  properties  about 
the  measurement  process  are  not  known,  but  statistical  properties  are  known  or  can  be  reasonably  well  modeled. 
Methods  using  PSDs  describe  the  performance  of  the  system  on  average,  and  any  optimization  of  a  system 
using  such  techniques  may  not  be  optimal  for  a  specific  instantiation  (i.e.,  a  matched  filter  would  be  best  for  an 
instantiation  where  the  object  and  any  interfering  signals  are  deterministically  known). 

In  remote  sensing  applications,  polarimeters  are  often  used  to  image  large  areas  of  similar  terrain,  and  they 
are  optimized  for  measurement  of  the  entire  data  set  as  opposed  to  one  particular  view.  The  fact  that  these 
imagers  measure  large  quantities  of  data,  and  have  prior  knowledge  of  signal  statistics  for  given  applications, 
makes  comparison  based  on  the  PSDs  a  reasonable  choice.  For  instance,  in  satellite  imaging  applications,  an 
existing  large  pool  of  data  allows  estimation  of  the  nominal  intensity  statistics  for  a  variety  of  wavelengths,  even 
if  the  polarization  PSDs  are  not  known. 

For  Stokes  polarimeters,  the  object  to  be  measured  can  be  represented  as  a  vector,  i.e.,  an  image  of  three 
(linear  polarimeters)  or  four  (complete  polarimeters)  Stokes  parameters.  The  power  spectrum  must  therefore  be 
represented  as  a  spectral  density  matrix.  The  spectral  density  matrix  (SDM)  is  defined  as  the  expected  value  of 
the  outer  product  of  the  Stokes  parameters  (S (/))  as  a  function  of  frequency, 

SDM(/)  =  <S(/)St(/)>  =  (||S(/)||2)  .  (4.1) 

The  SDM  contains  information  about  the  Fourier  transforms  of  the  auto-  and  cross-correlation  functions  among 
the  Stokes  parameters. 

The  goal  of  this  effort  is  to  develop  a  spectral  density  response  function,  analogous  to  the  transfer  function 
of  conventional  optics.  The  spectral  density  response  function  relates  the  SDM  of  the  object  Stokes  parameters 
to  the  spectral  densities  of  the  image  Stokes  parameters.  However,  the  measurement  process  for  modulated 
polarimeters  is  not  linear  shift  invariant  (LSI),  so  the  resulting  spectral  density  response  function  must  be 
able  to  characterize  the  magnitude  of  response  of  all  output  frequencies  for  a  given  input  frequency,  which  is 
somewhat  more  complicated  than  for  a  LSI  system.  As  an  example,  consider  a  pushbroom  polarimeters  such 
as  the  Multiangle  Spectro-Polarimetric  Imager  (MSPI)  J3HI  which  builds  up  an  image  one  line  at  a  time  as  the 


platform  moves  across  the  scene.  For  these  instruments,  spatial  variation  is  translated  into  the  temporal  response 
of  the  detector,  which  causes  spatial  variations  and  temporal  noise  to  combine  with  the  temporal  modulation 
strategy.  The  formalism  presented  here  can  be  applied  to  such  imagers. 

4.2  Spectral  density  response 

4.2.1  Derivation 

Consider  designing  an  imaging  polarimeter  for  an  application  where  the  Stokes  parameter  SDM  is  known.  In 
this  situation  it  is  desirable  to  predict  the  PSDs  of  the  output  Stokes  parameters  and  modify  the  polarimeter 
design  to  optimize  the  output  based  on  criteria  such  as  the  signal  to  noise  ratio  (SNR).  The  polarimeter  imaging 
operator  is  defined  as 

S[nx,  ny,  nt\  =  V  {S(:r,  y ,  t)}  ,  (4.2) 

where  S (x,y,t)  is  a  continuous  incident  set  of  object  Stokes  parameters  and  S[nx,ny,rit]  is  a  discrete  estimate 
of  the  incident  polarization  quantities.  The  polarimeter  operator  includes  both  the  physical  measurement  of  the 
incident  fields  and  the  reconstruction  from  that  measurement  into  the  final  estimation  of  the  scene.  This  operator 
is  useful  because  it  provides  a  description  of  the  entire  system  from  measurement  to  estimation  at  once,  instead 
of  first  choosing  a  measurement  scheme  and  then  the  estimation  algorithm.  Optimizing  the  polarization  operator 
for  an  imaging  task  will  potentially  lead  to  different  solutions  from  conventional  optimization  techniques.^ 

The  operator  V  is  developed  from  a  basic  description  of  the  system.  This  discussion  will  address  time 
modulated  polarimeters,  since  such  instruments  are  single  variate.  However,  the  approach  can  be  readily  applied 
to  more  dimensions,  such  as  spatially  modulated  polarimeters  (microgrids )P  or  a  system  modulated  in  both  space 
and  timep  and  is  directly  applicable  to  all  periodically  modulated  instruments.  Generalization  to  non-periodic 
polarimeter  modulation  schemes4' I  is  left  to  future  work.  The  optics  will  be  assumed  to  have  an  ideal  (flat) 
transfer  function  (infinite  aperture  diffraction  limited  optics),  but  the  detector  impulse  response  is  explicitly 
included  and  modeled  as  a  finite  integrator: 


h(t)  =  rect 


(4.3) 


where  rect(f)  =  1  for  |t|  <  1/2,  rect(f)  =  1/2  for  |i|  =  1/2,  and  0  otherwise.  The  sampling  function  is  defined  as 

W.(n,t)=6{t-nT),  (4.4) 

which  assumes  periodic  sampling  with  a  sampling  interval  of  T  and  sample  index  n.  The  flux  that  passes  through 
a  time- varying  analyzer  is  given  by 

P(*)=aT(W),  (4.5) 

where  a(f)  is  the  analyzer  vector  (the  first  row  of  the  Mueller  Matrix  of  the  polarimeter  opticsJ=Sl  This  flux  is 
then  sampled  by  the  detector 

/OO 

(h(t)  *  P(t))  IH(n,  f)df,  (4.6) 

-OO 

where  *  is  the  convolution  operator.  The  band-limited  data  reduction  algorithm  proposed  in9  is  described  by 


S[n]  =  w[n\  *  Z  i[n]as[n]/[n], 


(4.7) 


where  as[n]  is  a  sampled  version  of  the  analyzer  vector,  w[n]  is  the  reconstruction  window  (which  will  be  discussed 
in  detail  below),  and  the  analyzer  inversion  matrix  is 


^  1[n]  =  ( w[n }  *aj[n]as[n]) 


-1 


(4.8) 


Equations  4.6  4.7  and  4.8  constitute  the  full  mathematical  model  of  the  polarimeter  operator  using  band-limited 
reconstruction.  This  operator  can  now  be  used  to  compute  the  spectral  density  response.  This  analysis  will  be 
most  easily  accomplished  in  the  Fourier  domain  for  the  periodically  modulated  polarimeters  considered  here. 


First  we  note  that  the  integral  over  all  time  in  Eq.  |4.6|can  be  computed  by  evaluating  the  Fourier  transform 
of  the  integrand  at  /  =  0: 

(4.9) 


I[n}=  (h(f)P(fj)  *ni(n,/) 


/— o 


The  Fourier  transform  of  P(i)  is  found  by  taking  the  Fourier  transform  of  Eq.  4.5 

n/)  =  aT(/)*S(/). 


((Hf)  (aT(/)  *  S(/)))  *  III(?r,  /)) 


Substituting  Eqs.  T9  and  4.10  into  Eq.  |4.7|  yields 

S[n]  =  w[n ]  *  Z_1[n]as[n] 

A  choice  that  seems  reasonable  for  the  sampled  analyzer  vector  calculation  is 

/»00 

as[n\  =  /  (h(t)  *  a(t))Ul(n,t)dt, 


f= o 


(4.10) 


(4.11) 


(4.12) 


since  this  samples  the  analyzer  vector  the  same  way  the  flux  is  sampled.  Another  reasonable  choice  might  be 
sampling  the  analyzer  with  ideal  point  sampling.  The  point  of  this  discussion  is  that  this  parameter  is  free  for 
the  algorithm  designer  to  choose  to  best  fit  the  underlying  measurement.  In  fact,  the  sampled  analyzer  vector 
can  be  augmented  with  additional  signals  so  that  the  estimation  of  the  Stokes  parameters  is  orthogonal  to  certain 
known  sources  of  error;  the  authors  of®  augmented  their  conventional  data  reconstruction  matrix  (DRM)  with 
linear  gradients  to  force  the  system  estimation  to  ignore  slowly  varying  components  of  the  signal. 


Equation  4.11  becomes  more  manageable  when  we  make  some  simplifications  to  a (/)  and  S (/).  The  first 


simplification  is  to  assume  that  a (i)  has  a  period  B  and  can  be  written  in  terms  of  its  Fourier  serieJ^ 

oo  /  h  \  00 

a (t)  =  Y  cb  exp  (  2ni—t  J  =  Y  cb  exp  (2nifbt), 

b—  —  oo  ^  '  b——oo 

with  ft,  =  b/B.  The  Fourier  transform  of  this  Fourier  series  is  a  weighted  sum  of  delta  functions 

OO  /  t  \  OO 

a  (/)  =  Y  C6<V~r)=  cbS(f  -  fb). 


b— — oo 


(4.13) 


(4.14) 


b=— oo 


We  will  also  assume  periodicity  of  the  quantity  Z  1[n]a„[n],  since  we  have  already  assumed  periodicity  of  the 
analyzer  vector.  Likewise,  the  Fourier  series  of  this  function  is 


Z  1fnlas[nl 


=  Y  dp  exp  )  =  Y  dp  exp  (2niTfpn) 

p= 0  ^  p—0 


(4.15) 


where  P  is  the  period  of  Z_1[n]as[n],  p  is  an  integer,  and  fp  =  p/ PT.  Note  that  P  has  no  units  since  it  is  the 
period  of  a  discrete  set,  while  B  has  units  of  time  since  it  is  defined  in  the  continuous  domain.  The  summation 
only  has  P  elements,  i.e.,  the  number  of  measurements  taken  over  one  period  of  the  function  in  Eq.  |4.15| 
The  period  T  must  be  included  in  the  definition  of  fp  so  that  this  quantity  has  units  of  frequency.  Another 
simplification  made  here  is  the  introduction  of  the  quantity 


PU  -  fo)  =  S (/)<*(/  -  fo), 


(4.16) 


where  —  is  a  frequency-domain  function  that  represents  a  pure  sinusoidal  excitation  in  time  at  a  frequency 
given  by  f0  with  a  Fourier  domain  amplitude  described  by  S (/).  This  modification  is  made  so  that  the  response 
function  can  be  built  up  by  integrating  over  all  fQ  to  find  the  total  response  to  the  input. 


S(/)  =  P{f-  foWo. 


—  oo 


(4.17) 


The  function  (3(f  —  f0)  will  replace  S (/)  in  the  above  equations  for  the  derivation  of  the  spectral  density  response 
function. 


Using  the  simplifications  in  Eqs.  4.14  and  4.15  the  representation  of  the  sampled  flux  in  Eq.  |4.9| changes  to 


I[n]  =  Hf)  E  CW  ^  fb)  *  S(/)  *  III(n,  /) 


\b——oo 


(4.18) 


/= o 


The  convolution  of  S (/)  with  delta  functions  simply  shifts  the  argument  of  S(/)  in  frequency,  therby  creating 
side  bands  in  the  Fourier  domain— 


I[n]  =  Hf)  E  c?s(/  -  h)  *  ni(n,  /) 


\b——oo 


f= 0 


(4.19) 


When  examining  the  definition  of  I[n]  in  Eq. 
n.  The  Fourier  transform  of  III(n,  t)  over  t  is 


4.6 


we  find  that  only  the  sampling  function  HI  (n,i)  depends  on 


Tt^f  {III  (n,t)}  =  exp  (27 rinTf) , 

and  the  discrete  Fourier  transform  (DFT)  of  the  sampling  function  over  index  n  is 

N- 1 


1*  i-  /  ~  7  ”  \ 

Fn^k  {ni(n,/)|  =  E  exP  [Zmn  Tf  -  —  ) 

n— 0  '  *- 


1  —  exp  (27riNT f) 

1  -  exp  (2tt iT  [f  ~  j^]) 


(4.20) 


(4.21) 


Since  k  is  an  integer  and  does  not  have  units  of  frequency,  a  useful  variable  to  define  is  the  sampled  frequency 
fk- 


(4.22) 


which  relates  the  sampled  frequencies  back  to  the  excitation  frequencies.  Using  Eq.  |4.22|in  Eq.  |4.21|  yields 


in  (/*,/) 


1  —  exp  (2niNTf) 

1  -  exp(2? iiT  [f  -  fk])' 


(4.23) 


The  behavior  of  the  sampling  function  is  not  ideal;  the  Fourier  transform  of  the  ideal  sampling  function  would 
have  the  form  S(f  —  fk),  where  non-zero  reconstructed  values  would  only  occur  exactly  when  the  sampling 
frequency  equals  the  incident  frequency.  The  form  of  Eq.  |4.23|  allows  frequency  leakage  from  frequencies  where 
fk  7^  /•  Frequency  leakage  is  a  fundamental  sampling  issue,  but  is  specifically  discussed  here  to  fully  develop  a 
model  of  modulated  imaging  polarimetry.  Now  we  can  consider  that  the  sampling  delta  function  is  changed  to 
include  two  shifts  introduced  by  the  modulation  of  the  measuring  analyzer  vector  and  the  modulation  of  the  data 
processing  algorithm.  When  the  integers  b  and  p  are  zero,  the  sampling  occurs  exactly  at  fk  =  f  just  as  with 
a  conventional  imager.  As  b  and  p  change  they  allow  the  excitation  of  a  pure  sinusoidal  input  with  a  frequency 
/  to  create  a  response  in  the  system  at  a  frequency  given  by  fk  =  f  +  fb  +  fp ■  The  sampling  still  occurs  at 
specific,  discrete  intervals,  but  now  there  is  ambiguity  in  the  estimation  because  multiple  excitations  can  create 
overlapping  responses. 

There  are  two  aspects  to  consider  when  dealing  with  sampling  artifacts  in  polarimeters:  aliasing  and  frequency 
leakage.  For  the  examples  used  in  this  discussion,  the  effects  of  aliasing  will  be  shown  in  the  following  section. 
The  response  of  the  real  system  to  frequencies  where  f  ^  fk  is  given  by  Eq.  |4.23|  which,  for  /  sufficiently 
close  to  fk,  can  be  approximated  with  the  function  Ul(fk,f)  =  sine (T(f  —  fk)),  where  the  sine  function  allows 
leakage  from  a  band  in  /  to  a  particular  value  fk-  The  behavior  when  /  =  fk  is  still  ideal,  but  for  frequencies 
fk- 1  <  /  <  fk  there  is  leakage.  The  frequency  leakage  behavior  can  be  included  by  convolving  the  ideal  results 
with  the  sine  function  that  models  the  leakage,  and  will  be  ignored  for  now. 


Eq. 


With  the  Fourier  transform  over  /  complete,  now  we  need  to  take  the  discrete  Fourier  transform  (DFT)  of 
over  the  index  n.  First,  substituting  the  Fourier  series  of  the  sampled  analyzer  vector  and  Z  yields 


4.11 


p- 1 


S  W  =  w[n]  *  E  dp  exp  f  -  j  I[n]. 

p—0  '  ' 

Taking  the  DFT  of  S[n]  and  using  properties  of  delta  function^H  provides 

p-i 

S [fk]  =  E  dPw[fk]I\fk  ~  fP]- 


(4.24) 


(4.25) 


p—0 


Now  that  the  total  polarimeter  operator  is  represented  in  the  Fourier  domain,  the  spectral  density  response 
can  be  calculated.  The  set  of  power  spectra  for  the  estimated  Stokes  parameters  (or  spectral  density  response) 
is  defined  as  the  Fourier  transform  of  the  auto  correlation  function.  This  can  be  written  as 

Wk]  =  {s  [n]  *  S  [n]}  =  kfk]&[fk],  (4-26) 

where  t  represents  the  Hermitian  adjoint  and  *  is  the  correlation  operator.^  This  quantity  can  be  calculated  in 
the  Fourier  domain  as 


E  [h\ 


P—1  oo 

E  E  Hf-fk+fP)h*(f-fk+fP')\\ti[fk]\\2 

p,p'—0  b,b'=— oo 


(4.27) 


where  fx  =  fk-fP+fb,  fi  =  fk~fP’+fb',  and  the  object  SDM  =  (s (/  -  /i)S(/  -  /2)t^)  describes 

the  cross  spectral  density  of  the  Stokes  parameters  and  the  cross  spectral  density  of  differently  modulated  Stokes 
parameters.  Note  that  *  denotes  the  non  transposed  complex  conjugate  of  a  vector.  Equation  |4.27|  provides  a 
solution  to  the  problem  that  we  set  out  to  solve  originally:  given  a  known  SDM  for  the  object  Stokes  parameters, 
what  is  the  spectral  density  response  of  the  estimated  Stokes  parameters?  Using  these  methods,  and  with  a 
known  estimate  of  the  incident  Stokes  parameters,  the  spectral  density  response  R[/fe]  can  be  directly  computed. 


4.2.2  Illustration 


To  illustrate  the  concepts  leading  up  to  Eq.  4.27  consider  the  case  of  a  rotating  analyzer  polarimeter,  which  can 
be  described  by 

1 


a(i)  =  \ 


C0S(27 Tyt) 
sin(27r^t) 
0 


(4.28) 


where  fs  is  the  sampling  rate  of  the  sensor.  The  Fourier  transform  of  the  rotating  analyzer  vector  is  shown  in 
Fig.  |4.1|  We  shall  also  assume  that  the  system  performs  ideal  periodic  point  sampling,  such  that  h(t)  =  5(t),  or 


Kf)  =  i. 


Even  though  the  aim  of  this  paper  is  to  describe  the  operation  of  a  modulated  polarimeter  on  stochastic 
signals,  the  steps  from  the  previous  subsection  are  more  easily  visualized  using  a  deterministic  signal.  In  this 
illustration,  we  assume  that 

r  V2 1 


S  (t) 


1 

1 


sinc2(f), 


(4.29) 
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Figure  4.1:  Analyzer  vector  for  a  rotating  polarizer  polarimeter  in  the  Fourier  domain.  Since  the  analyzer  is  comprised 
of  periodic  functions,  in  the  Fourier  domain  it  is  simply  comprised  of  shifted  delta  functions. 


Figure  4.2:  The  components  of  the  Stokes  parameters  so  -  So  modulate  the  carriers  a o  -  02  to  produce  the  final  sampled 
signal  I[n].  The  demodulation  process  is  composed  of  remodulation  using  ao  -  a 2,  then  unmixing  and  equalizing  the 
channels  using  Z  .  The  colors  allow  us  to  track  the  contributions  from  s 0  -  S2  through  the  modulation  and  demodulation 
processes. 


The  steps  in  Eq.  |4.5|  -  Eq.  |4.8|  in  the  time  domain  are  shown  in  Fig.  |4.2|  and  the  demodulation  steps 
are  depicted  in  the  frequency  domain  in  Fig.  |4.3|  with  an  arbitrary  low-pass  filter.  The  figure  shows  how  the 
information  from  each  Stokes  parameter  channel  maps  into  the  components  of  the  final  signal. 

The  computation  of  the  resulting  spectral  densities  for  s0>  si>  and  s2  given  in  Eq.  |4.27|is  shown  in  Fig.  |4.4| 
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Figure  4.3:  In  the  frequency  domain,  the  three  baseband  signals  corresponding  to  the  Stokes  parameters  are  placed  into 
orthogonal  channels  by  the  carriers.  In  the  demodulation  process,  each  of  the  carriers  is  used  to  place  the  corresponding 
desired  side  band  at  baseband.  A  filter  is  used  to  isolate  the  desired  signal  to  create  the  final  estimate. 


Figure  4.4:  Graphical  representation  of  the  estimation  of  the  spectral  density  response  function  before  filtering,  showing 
the  contributions  due  to  all  three  parameters.  The  total  spectral  density  for  each  parameter  is  obtained  by  summing  all 
contributions.  Note  that  the  signal  for  S3  is  shown  with  a  dashed  red  line,  (a)  Reconstruction  of  so-  (b)  Reconstruction 
of  si.  (c)  Reconstruction  of  S2. 


f/f.  S/S,  S/S, 


Figure  4.5:  The  assumed  object  SDM  for  this  exercise.  The  diagonal  elements  correspond  to  the  spectral  density  for  (a) 
So,  (e)  Si,  and  (i)  S2.  The  off  diagonal  elements  represent  the  cross-spectral  densities,  for  example  (b)  is  the  cross-spectral 
density  between  Si  and  So- 


4.3  Example:  Computation  of  the  Wiener  filter  for  modulated  polarimeters 


Finally,  we  assume  that  the  SDM  for  the  object  to  be  measured  is  band-limited  and  described  according  to  the 
equation 
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■  1 
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(4.31) 


The  object  SDM  is  assumed  to  be  uncorrelated  for  differently  modulated  Stokes  parameters  /1  7^  and  the 

4.5  Note  that  because  03(1)  =  0,  the  S3  information  is  lost,  thereby  allowing  us  to 
3x3  matrices.  It  would  also  be  possible  to  retain  the  dimensionality,  but  because  Z 


matrix  is  shown  in  Fig. 
reduce  Z_1  and  R(/)  to 


would  no  longer  be  square,  we  would  have  to  resort  to  using  a  pseudoinverse  instead. 


The  classical  DRM  method28  (referred  to  here  as  “the  sliding  pseudo  inverse”)  and  the  ideal  band-limited 
filters-^  have  been  presented  in  the  literature,  but  for  practical  applications  neither  filter  performs  adequately 
for  imaging  polarimetersP^  The  sliding  pseudo  inverse  produces  strong  artifacts  at  edges  in  the  object,  while 
the  ideal  band-limited  filter  only  performs  well  on  band-limited,  noiseless  data;  estimation  from  noisy  or  non- 
band-limited  data  often  produces  strong  ringing  artifacts  in  the  space-time  domain.  A  more  “visually  pleasing” 
band-limited  filter  would  be  one  that  has  a  smooth  transition  to  zero,  rather  than  the  sharp  cutoff  that  the  ideal 
band-limited  filter  has.  One  approach  to  building  a  smooth  band-limited  filter  is  the  Wiener  filter.  The  Wiener 
filter  is  well  studied;  one  in-depth  version  of  the  Wiener  filter  and  Wiener-Helstrom  theory  is  in.°°  This  filter  is 


derived  for  non  correlated  white  noise  using  the  metric  of  least  square  error.  The  Wiener  filter  is  used  here  only 
as  an  example  to  illustrate  the  application  of  the  current  theory,  since  it  is  a  well-defined,  optimal  linear  filter. 
Other  filters  or  nonlinear  reconstruction  methods  may  be  more  desirable  for  a  particular  application,  and  the 
spectral  density  response  formalism  developed  above  can  be  used  for  those  as  well. 

The  Wiener  filter  is  designed  to  solve  the  estimation  problem  described  by 


y(t)  =  w(t)*g(t), 


(4.32) 


where  g(t)  is  the  measured  signal 


g(t)  =  s(t)  +  n(t). 


(4.33) 


In  Eq.  4.33  s(t)  is  the  true  signal  and  n(t)  is  additive,  uncorrelated  noise.  The  Wiener  filter  is  found  by  solving 

w(t)  ( g{t )  *  g(t))  =  s(t )  *  g(t),  (4.34) 

and  in  the  Fourier  domain  the  Wiener  filter  is 

IK/)II2 


K/)  = 


K/)ii2+iH/)ir 


(4.35) 


The  Wiener-Helstrom  filter,  which  makes  use  of  knowledge  of  the  transfer  function  h(f)  of  the  system  (also 
derived  ino(3),  can  be  calculated  using 


™(f) 


h{f ) 


ls(/)l|2  +  lln(/)l|2 


(4.36) 


The  estimator  in  Eq.  |4.11|  unfortunately  has  two  instances  of  the  filter,  the  second  of  which  is  contained  in 
the  calculation  of  Z.  The  analytic  calculation  of  a  filter  that  minimizes  square  error  is  difficult  with  the  filter 
included  in  the  calculation  of  the  inverse.  However,  if  we  impose  some  loose  constraints  on  ui(f),  the  problem 
simplifies,  and  we  can  ignore  the  fact  that  ui(f)  is  included  in  the  calculation  of  Z.  Consider  the  portion  of  the 
operator  inside  of  the  inverse: 

g-\t)  =  ( w(t )  *  A(t)AT(f))  1  .  (4.37) 

To  ignore  the  w(t)  in  the  definition  of  Z(t),  the  condition 

/OO 

(A(f')A1  (t1))  dt1  (4.38) 

-OO 

must  be  satisfied.  If  we  assume  that  A(f)  is  periodic,  it  can  be  rewritten  in  terms  of  its  Fourier  series 


aj(t)  =  ^  zjm  exp  . 

n=— oo 

Using  Eq.  |4.39|in  the  right  side  of  Eq.  |4 . 38 1  yields 


OO  OO 

Zjj>  =  ^2  ^2  ZjmZj'm’Slm  +  m']. 


(4.39) 


(4.40) 


m' ——oo  m=— oo 


Using  Eq.  |4.39|in  the  left  side  of  Eq.  |4.38|  yields 

OO  OO  / 

)  =  ^2  X]  zimZj'm'W  (: t )  *  exp  f  27 u 

m'=— oo  m=— oo  ^ 
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(4.41) 


Taking  the  Fourier  transform  of  Eq.  |4.41|  yields 

OO  OO 

Zjj'(f)=  zjmZj'm'W{f)5[  f  - 


m  — — oo  m——oo 


m  +  m! 


(4.42) 


For  Eqs.  4.40  and  4.42  to  be  equivalent,  the  restrictions  imposed  on  the  filter  are 


1,  for  /  =  0 

w{f)=  {  0,  for  /  = 

arbitrary,  otherwise 


(4.43) 


These  restrictions  allow  us  to  consider  a  large  set  of  inverse  operators  for  which  the  calculation  of  Z  1  yields 
identical  results,  and  we  can  ignore  w(f)  in  its  calculation. 

For  the  moment,  assume  the  impulse  response  functions  associated  with  the  system  optics  and  detector  are 


ideal.  With  the  conditions  imposed  on  w(f)  in  Eq.  4.43  as  well  as  ideal  point  impulse  response  functions,  the 
unfiltered  estimate  of  the  Stokes  parameters  in  the  Fourier  domain  can  be  written  as 


P—1  OO 

sifki  =  Y  Y  dPcbS(f- h), 

p—0  b=— oo 


(4.44) 


with  fi  and  f2  defined  in  Section  4.2.1  Using  this  g [/*.]  in  the  calculation  of  the  Wiener  filter  without  noise 
provides 

J2P,bzp,icb  (s(/- 


Wi[fk\  =  — 

Ep.p',6,6'  ZP,iCI  (S(/  -  /l)S*(/  -  /2))  c *b,Zfp/  i 

Rewriting  this  in  terms  of  the  cross-spectral  density  matrix  defined  in  Eq.  |4.27| 


i>i[fk]  = 


E„E6sViE*(/;/-/i) 


E p, p',6,6'  zP,icJMf  -  fi;  f  -  h K'4',i 


(4.45) 


(4.46) 


Under  the  assumption  of  completely  uncorrelated  Stokes  parameters,  the  Wiener  filter  calculation  simplifies  to 
the  form 


Wi[fk\  = 


Si(f ) 


EP,  P',6,6'  zP,icjR(f  ~  A;  f  ~  h  )cl'zl',i 


(4.47) 


Assuming  that  the  Stokes  parameters  are  completely  uncorrelated  with  one  another  is  a  poor  assumption  from 
a  physical  perspective.  Even  if  the  polarization  channels  of  the  Stokes  parameters  are  uncorrelated  with  each 
other,  So  must  have  some  degree  of  correlation  with  Si,  S2,  and  S3.  This  is  because  physically  realizable  signals 
always  have  DOLP  less  than  1.  However,  this  assumption  led  us  to  a  form  of  the  Wiener  filter  that  compared 
well  analytically  to  the  known  result  for  conventional  imaging,  and  is  useful  to  discuss  this  result  further. 

Consider  the  system  to  be  corrupted  by  additive  white  noise,  changing  the  unfiltered  estimation  to  be 


P—1  OO 

g  n[fk\  = 

p— 0  b=— 00 


<£s(/- A) +  «(/-/*  +  /„) 


(4.48) 


The  spectral  density  response  found  with  only  noise  (zero  contribution  from  Stokes  parameters)  will  from  now 
on  be  referred  to  as  the  noise  response  function  (NRF), 


dP«(/  ~  fk  +  ~  fk  +  /P')dp'  |  . 

1  p,p' 


NRF(/)  =  diag 


(4.49) 


where  the  quantity  n(f  —  fk  +  fP)n*(f  —  fk  +  fp '),  describing  the  power  spectrum  of  the  noise,  must  be  known 
or  estimated  for  the  application.  This  definition  is  similar  to  the  required  knowledge  of  the  SDM  for  the  Stokes 
parameters.  With  added  noise,  the  Wiener  filter  calculation  becomes 


Wi[fk] 


Si(f) 


Ep, P',6, 6'  zP,iclMf  —  /i ;  /  —  h)c*b'Zp>  i  +  NRFi(f) 


(4.50) 


Once  more  including  the  transfer  functions  of  the  detector,  the  unfiltered  spectral  density  response  is  described 
by 


l|g[/fc]||2  =  £  £M/  -  fk  +  fP)h*(f  -fk  +  fp') 

p,pf  b,b' 

x  dpcjE(/-/i;/-/2)Cb,d*,, 


(4.51) 


and  the  Wiener-Helstrom  filter  is 


i[fk]  = 


h*(f)  Si(f) 


II  9i[fk 


NRFi(f) 


(4.52) 


The  Wiener-Helstrom  filter  has  the  advantage  that  it  can  attempt  to  recover  information  that  was  deterministi¬ 
cally  attenuated  due  to  the  system  transfer  function. 


4.4  Discussion 

4.4.1  Infinitesimally  small  detector  integration  time 


Spectral  Density  Response  When  the  statistical  object  is  measured  by  the  rotating  analyzer  system  and  subse¬ 
quently  estimated  as  described  above,  the  estimation  results  in  a  weighted  superposition  of  shifted  copies  of  the 
object  SDM  (Eq.  4.27 1.  Figure  4.4  shows  that  the  estimation  of  each  Stokes  parameter  is  correctly  centered  at 
the  low  frequency,  but  there  are  also  higher  frequency  copies  of  all  of  the  other  Stokes  parameters  that  interfere 
with  the  estimation;  these  copies  must  be  filtered  out  with  a  low  pass  filter  w,  which  will  replace  the  filter 
with  uniform  frequency  response.  The  spectral  density  response  function  using  a  uniform  response  filter  allows 
us  to  find  the  Wiener  filter  that  is  optimal  for  estimating  objects  that  have  the  same  statistical  properties  as 
the  object  we  have  modeled.  The  Wiener  filter  calculation  is  shown  in  Eq.  |4.47|  graphically  the  Wiener  filter 
calculation  for  the  *th  Stokes  parameter  will  be  performed  by  dividing  Ru(f )  (the  diagonal  elements  in  Fig.  4.5  I 


by  the  corresponding  spectral  density  response  obtained  from  Fig. 


The  resulting  Wiener  filters  are  shown 


in  Fig.  |4.6|  These  filters  are  band-limited  filters.  In  fact,  the  Wiener  filter  for  a  system  that  has  no  noise  and  no 
overlapping  signals  converges  to  the  ideal  rectangular  band-limited  filter.  However,  the  Wiener  filter  smoothly 
approaches  zero  in  the  presence  of  overlapping  signals,  which  allows  the  filter  to  minimize  undesirable  ringing 
artifacts  in  the  temporal  domain  in  an  objectively  determined  way. 


We  can  also  calculate  a  Wiener  filter  for  functions  that  are  assumed  to  be  completly  correlated.  We  will 
assume  that  the  Stokes  parameters  are  described  by  the  same  power  spectra  and  that  differently  modulated 
Stokes  parameters  are  uncorrelated,  but  that  each  of  the  Stokes  parameters  are  completly  correlated  with  each 
other.  The  graphical  representation  of  the  spectral  density  response  function  estimation  for  this  case  is  shown 
in  Fig.  |4.7[  The  effect  of  a  given  cross-term  is  cancelled  out  by  the  effect  of  its  complex  conjugate,  which  is 
also  added  to  the  spectral  density  response  estimation.  This  is  represented  by  the  magenta  line  at  the  bottom 
of  the  figure.  Thus,  as  long  as  differently  modulated  parameters  remain  uncorrelated,  the  estimation  of  the 
spectral  density  response  remains  the  same,  whether  the  Stokes  parameters  are  correlated  or  uncorrelated  with 
each  other. 


Frequency-domain  map  of  the  SDM  To  further  our  understanding  of  the  a  system’s  spectral  response,  we 
note  that  the  SDM  in  Eq.  |4.27|  can  be  represented  as  a  frequency-domain  map  relating  the  spectral  density 
response  of  the  estimated  Stokes  parameters  at  different  sampled  frequencies  to  the  estimated  spectral  density 
response  of  the  continuous  Stokes  parameters.  Figure  |48|  is  an  example  of  elements  i?ocb  fin,  and  R22  of  the 
estimated  unfiltered  SDM  for  the  object  represented  in  Eq.  |4.31[  In  this  case,  unfiltered  means  that  the  filter  in 
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Figure  4.6:  Wiener  filters  calculated  for  the  assumed  Stokes  parameters  power  spectra  being  measured  by  the  rotating 
analyzer  polarimeter.  The  Wiener  filters  follow  the  signal  to  noise  ratio  of  the  estimation  as  a  function  of  frequency. 
These  filters  are  band-limited,  but  approach  the  cutoff  smoothly  to  avoid  undesirable  temporal  ringing  artifacts. 


Figure  4.7:  Graphical  representation  of  the  spectral  density  response  function,  assuming  a  filter  with  uniform  frequency 
response,  for  (a)  So,  (b)  Si,  and  (c)  S2,  and  completely  correlated  Stokes  parameters.  Curves  in  blue,  green,  and  red  are 
components  originally  due  to  So,  Si,  and  S2,  respectively.  The  effect  of  the  cross  terms  is  represented  by  the  magenta 
line  at  the  bottom  of  the  figure.  The  total  spectral  density  for  each  parameter  is  obtained  by  summing  all  contributions. 


the  reconstruction  is  taken  as  w  =  1  over  the  whole  frequency  domain.  Since  we  assumed  that  different  Stokes 
parameters  are  uncorrelated,  elements  Rij  are  zero  for  i  7^  j.  In  these  figures,  the  vertical  axis  is  the  frequency  of 
the  continuous  signal,  whereas  the  horizontal  axis  is  the  sampled  frequency.  The  total  spectral  density  response 
obtained  with  the  superposition  of  shifted  elements  from  the  object  SDM  represented  in  Fig.  |4.4|  is  a  cross 
section  of  Fig.  |4.8|  along  the  line  /fc  =  0.  Therefore,  the  SDM  frequency-domain  map  is  a  more  general  way  of 
representing  the  system’s  response,  since  it  contains  the  response  at  different  sampling  frequencies. 


The  secondary  maxima  in  Fig.  4.8  constitute  artifacts  that  limit  our  ability  to  reconstruct  the  SDM  of  the 
input  signal.  However,  the  Wiener  filter  presented  in  Section  |4.3|  allows  us  to  filter  out  the  higher  frequency 
replicas  in  the  unfiltered  signal  and  can  also  be  represented  as  a  frequency-domain  map.  Figure  |4.8|  is  the 
modulus  square  of  the  Wiener  filter  frequency  map  for  each  Stokes  parameter.  Just  as  for  the  SDM  frequency 
maps,  a  cross  section  of  this  figure  at  fk  =  0  results  in  the  Wiener  filter  results  in  Fig.  |4.6|  Applying  this  filter 
during  the  estimation  of  the  Stokes  parameters  represented  in  Eq.  |4.27|  results  in  the  filtered  signal  that  yields 
the  SDM  in  Fig.  4.8  The  filtered  SDM  elements  closely  resemble  a  tri z(f)  along  the  diagonal,  which  would  be 
the  spectral  density  response  function  of  an  ideal  system  with  no  polarization  artifacts. 


4.4.2  Numerical  comparison  between  filter  types 

In  this  section  we  quantitatively  compare  filter  performance  for  the  uncorrelated  example  using  the  mean  squared 
error  metric  (under  which  the  Wiener  filter  is  derived).  The  filters  we  use  for  this  section  are  the  sliding  pseudo 
inverse,  the  ideal  band-limited  filter,  and  the  polarimetric  Wiener  filter.  The  frequency  response  functions  of  the 
filters  are  shown  alongside  the  So  spectral  density  response  in  Fig.  |4.9|  To  begin  the  numerical  comparison,  we 
generated  an  object  with  the  statistical  properties  that  were  assumed  in  Eq.  |4.31|  To  do  this,  many  instantiations 
of  white  noise  were  added  to  the  Stokes  parameters  of  the  object  to  build  the  desired  SDM  with  a  signal-to-noise 


Figure  4.8:  Frequency-domain  map  for  the  estimated  (a)  Roo,  (b)  Ru,  and  (c)  R22  unfiltered  elements  of  the  SDM. 


Figure  4.9:  An  illustration  of  the  frequency  repose  of  (a)  the  sliding  pseudo  inverse,  (b)  the  ideal  band-limited  filter,  and 
(c)  a  Wiener  filter  (smooth  band-limited  filter)  compared  to  the  graphical  representation  of  the  spectral  density  response 
of  the  So  channel. 


ratio  of  20  dB.  A  single  instantiation  of  the  white  noise  is  shown  in  Fig.  4.10  Note  that  this  simulation  does  not 
guarantee  physical  results  (degree  of  polarization  less  than  1);  physicality  could  be  enforced  by  adding  a  constant 
offset  to  the  So  object.  However,  the  purpose  of  this  discussion  is  to  describe  the  method  of  application  of  these 
techniques. 

The  matrix  obtained  after  1000  iterations  is  shown  in  Fig.  |4.11|  which  is  similar  to  the  ideal  SDM  in  Fig.  |4.5 
but  with  added  noise.  This  SDM  approximates  the  desired  object,  as  the  off  diagonal  elements  are  approximately 


within  2.5%  of  zero.  The  object  is  measured  and  estimated  by  the  polarimeter  using  the  various  filters  in  Fig.  4.9 
The  mean  square  error  for  the  estimation  across  all  instantiations  is  shown  in  Fig.  |??[  The  two  bottom  rows 
are  quite  similar  and  show  an  overall  better  performance  than  the  ideal  band-limited  filter.  The  Wiener  filter 
distribution  is  slightly  flatter  than  that  of  the  sliding  pseudo  inverse,  and  continues  decreasing  until  the  edge  of 
the  frequency  region  covered  by  the  original  Stokes  parameters  power  spectra,  whereas  the  sliding  pseudo  inverse 
reaches  a  minimum  and  then  rises  at  the  edge  of  the  considered  frequency  region.  Hence,  the  Wiener  filter 
outperforms  either  of  the  other  two  filters,  which  was  expected  since  it  was  designed  for  optimal  performance 
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Figure  4.10:  Single  instantiation  of  the  white  noise  for  (a)  So,  (b)  Si,  and  (c)  S2 ■  The  noise  is  not  guaranteed  to  produce 
a  physically  realizable  signal  since  the  object  SDM  employed  in  this  problem  does  not  enforce  that. 
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Figure  4.11:  Calculated  SDM  for  the  simulated  Stokes  parameters  with  added  white  noise.  Here  the  SDM  (diagonal) 
elements  for  (a)  So,  (e)  Si,  and  (i)  S2  are  approximately  the  desired  triangle  squared  functions.  The  off  diagonal  cross 
spectral  density  elements  ideally  should  be  zero,  and  they  are  within  2.5%  of  the  correct  answer. 


under  the  MSE  metric.  Note  that  this  analysis  is  not  meant  to  demonstrate  the  superiority  of  the  Wiener 
filter,  rather  to  demonstrate  that  the  theory  presented  here  results  in  a  Wiener  filter  that  is  consistent  with 
expectations. 


4.4.3  Finite  detector  integration  time 

Now  consider  a  case  where  there  is  a  non-unity  detector  transfer  function,  such  that  the  detector  integration 
time  has  a  duty  cycle  of  1  and  is  described  in  time  as 


hit)  =  rect  (  - 


In  the  Fourier  domain,  the  detector  transfer  function  is  written  as 


(4.53) 


h(f)  =  sine 


/, 


(4.54) 


The  first  zero  of  the  detector  transfer  function  for  this  problem  is  then  at  fs,  while  we  are  concerning  ourselves 
with  signals  that  oscillate  at  rates  of  less  than  0.5 f3.  Figure  |4.12|  shows  the  graphical  representation  of  the 
spectral  density  response  i?oo  measured  with  a  detector  with  finite  integration  time,  together  with  the  detector 
transfer  function,  compared  to  the  ideally  measured  spectral  density  response  and  the  difference  between  these 
two  responses.  We  observe  that  the  So  channel  has  been  symmetrically  affected  by  the  detector  transfer  function 
since  they  have  the  same  center  of  symmetry,  but  the  high  frequency  shifted  Si  and  S2  channels  have  not  been 
symmetrically  affected.  Many  conventional  processing  techniques  for  modulated  polarimetry  implicitly  ignore 
this  non  symmetrical  effect,  since  conventional  calibration  is  done  at  the  zero  frequency.  The  lower  the  duty  cycle 
of  the  sampling  detector,  the  less  effect  this  asymmetry  has  since  the  detector  is  closer  to  ideal  point  sampling. 
Had  the  modulation  scheme  been  at  a  higher  frequency,  the  shifted  copies  would  be  further  away  from  the 
symmetrical  point  and  the  estimation  would  be  further  degraded.  The  Wiener-Helstrom  filter  can  be  computed 
in  this  case,  although  for  this  particular  polarimeter  little  advantage  is  gained.  The  effect  of  the  detector  transfer 
function  being  included  in  the  simulation  had  a  maximum  effect  of  about  2.5%  at  certain  frequencies,  and  this 
translates  to  having  a  maximum  effect  on  the  filter  of  less  than  0.5%  at  those  same  frequencies.  The  Wiener- 
Helstrom  filter  is  shown  in  Fig.  |4.12  as  well  as  the  difference  between  the  Wiener  and  the  Wiener-Helstrom 


Figure  4.12:  (a)  Graphical  representation  of  the  spectral  density  response  of  a  system  with  a  detector  with  finite  integration 
time  of  duty  cycle  1.  The  modulus  square  of  the  detector  transfer  function  is  shown  for  reference,  (b)  Graphical 
representation  of  the  spectral  density  response  using  an  ideal  detector  with  uniform  frequency  response,  (c)  The  difference 
between  the  spectral  density  response  with  an  ideal  detector  and  the  detector  with  finite  integration  time  in  (a),  showing 
the  non  symmetrical  effect  on  the  high  frequency  copies. 


filters.  This  difference  is  quite  small,  but  for  other  applications  there  may  be  differing  results.  The  importance 
of  the  detector  transfer  function  will  increase  for  applications  where  the  modulation  scheme  introduces  several 
side  bands,  such  as  the  MSPI  polarimeter  described  by  Diner,  et  Here,  the  advantage  of  calculating  the 
Wiener-Helstrom  filter  may  lie  in  helping  to  control  channel  cross  introduced  by  the  spatial  bandwidth  (converted 
to  temporal  bandwidth  by  the  moving  platform)  of  the  scene  being  measured. 

The  techniques  presented  here  provide  a  way  to  study  the  effect  measurement  has  given  models  of  the  SDM 
of  objects  to  be  measured.  The  development  of  the  optimal  Wiener  filter  for  a  particular  application  is  highly 
dependent  on  the  accuracy  of  the  object  models.  However,  even  without  accurate  object  models,  rough  guesses 
of  the  object  model  may  yield  effective  first-guess  filters  than  can  be  iteratively  improved  upon  as  empirical 
data  become  available.  Intuitively  we  know  that  a  smooth  band-limited  filter  should  provide  better  results  for 
periodic  analyzer  vectors  than  the  ideal  band-limited  or  sliding  pseudo  inverse:  it  is  band-limited  to  prevent 
channel  crosstalk  and  smooth  to  prevent  ringing  in  the  space  time  domain.  The  optimal  Wiener  filter  will 
perform  better  than  the  other  linear  filters  in  terms  of  the  mean  squared  error  for  a  given  well  defined  set  of 
object  statistics,  because  the  Wiener  filter  is  derived  by  minimizing  this  metric.  The  detection  task  drives  the 
selection  and  optimization  of  the  filter.  Therefore,  for  any  given  application  other  filters  than  the  Wiener  filter 
may  be  more  appropriate.  The  approach  developed  here  allows  for  the  comparison  of  the  performance  of  different 
polarimeter  modulation  schemes  and  processing  algorithms  simultaneously  given  a  defined  application  (object 
SDM).  The  spectral  density  response  presented  can  be  employed  to  perform  objective  comparisons  between 
different  polarimeter  designs  for  suitability  for  an  imaging  task. 

5.  CHANNELED  POLARIMETERS 

In  the  class  of  channeled  or  modulated  polarimeters,  the  polarization  states  that  define  the  measurement  are 
modulated  either  spatially,  temporally,  or  spectrally P  Every  one  of  those  harmonic  modulations  will  split  the 
information  in  the  corresponding  Fourier  domains,  creating  weighted  copies  of  the  Fourier  Transform  of  the  data 
at  the  modulating  frequencies.  These  multiplexed  copies  are  called  channels.  Oka  and  his  coworkerdMI51  °2  have 
popularized  the  design  concepts  that  go  into  making  a  channeled  system,  which  were  then  further  developed  by 
Hagen, ^  KudenovEMHl  and  their  coworkers  before  going  into  wide  use  by  many  others.  In  this  paper  we  will 
introduce  a  toolkit  to  describe,  analyze  and  optimize  such  systems,  and  we  will  use  it  to  investigate  channeled 
polarimeters  from  the  literature  to  show  how  they  can  be  improved. 

There  has  been  a  number  of  proposed  channeled  systems  in  the  past*40  45  sl  °4  whose  designs  and  the  cor¬ 
responding  reconstruction  techniques  were  derived  by  hand.  By  limiting  the  number  of  parameters  available 
during  design,  the  likelihood  of  a  suboptimal  design  is  inadvertently  increased.  Lemaillet  et  al°L°  proposed  a 
way  to  optimize  a  spectrally  channeled  system  by  introducing  linear  algebra  methods  to  map  the  information. 
Their  effort,  however,  focused  on  the  particulars  of  one  kind  of  system  and  stopped  short  of  providing  a  complete 
solution  to  deal  with  any  channeled  polarimeter.  This  paper  describes  the  generalized  methods  that  can  be  used 
to  model  channeled  information  mapping  and  guide  the  reconstruction. 

A  great  advantage  posed  by  a  channeled  system  is  the  possibility  of  constructing  a  snapshot  system  that  is 
able  to  simultaneously  record  information  pertaining  to  different  polarization  information  channels.  This  removes 
the  need  for  complex  image  registration  that  would  be  required  in  a  temporally  modulated  system.  In  terms 
of  object  bandwidth,  a  snapshot  channeled  system  favors  temporal  resolution  at  the  cost  of  introducing  stricter 
band  limit  constraints  in  other  domains. 

We  consider  a  common  2D  FPA  detector,  which  enables  access  to  up  to  two  modulation  types  to  be  mapped 
onto  the  two  orthogonal  axes.  In  addition  to  having  no  modulation,  we  can  also  map  spatial  and  spectral 
modulations  into  either  x-  or  y- axes  of  the  detector.  Although  the  methods  introduced  here  are  general  enough 
to  be  used  with  any  channel  structure  on  any  orthogonal  coordinates  systems,  this  paper  will  focus  on  Cartesian 
coordinates,  implying  that  the  channels  lie  on  a  rectangular  grid. 

For  the  sake  of  completeness,  we  will  also  consider  temporal  modulation.  Such  a  system  will  obviously 
lose  its  snapshot  nature,  but  it  is  conceivable  that  some  middle  ground  solution  could  be  found,  whereby  a 
very  limited  number  of  temporal  measurements  are  made  with  intent  of  balancing  the  resolution  loss  among 
all  possible  modulation  dimensions J3HI  Thus,  using  a  conventional  detector  will  allow  us  to  split  polarization 


information  into  a  three  dimensional  structure  of  channels  that  we  can  manipulate  to  reconstruct  the  polarization 
information.  It  is  useful  to  recognize  that  when  we  have  a  small  number  of  temporal  measurements,  the  resultant 
temporal  frequency  channels  may  contain  as  few  as  one  data  point.  In  those  cases  it  may  be  more  prohibitive 
and  unnecessary  to  work  with  the  data  in  the  Fourier  domain.  We  can  instead  use  measurements  as  information 
“channels”  themselves  with  a  clear  benefit  that  they  will  contain  modulation  information  more  compactly. 


5.1  Sinusoidal  Channel  Splitting 

Every  sinusoidal  modulation  in  either  ai  or  gj  splits  the  element  information  in  m,7-  into  two  channels  at  certain 
frequencies  within  the  Fourier  domain  of  the  modulation.  For  the  available  modulation  dimensions  of  x,  y,  a 
(wavenumber)  and  t.  we  will  call  the  corresponding  frequency  dimensions  £,  77,  r  (OPD)  and  v.  We  will  only 
show  the  relevant  equations  for  the  pair,  since  all  others  can  be  obtained  by  taking  any  other  sequential 
pairing  of  the  above  mentioned  dimensions.  The  following  transform  pairs  are  well  knownP^ 


l(z)  <- 

(5.1a) 

cos(27r^jai)  <- 

+  &)  +  <5(£  —  &)], 

(5.1b) 

sin(27r£ja;)  <- 

->•  *[««  +  &) -*«-&)]• 

(5.1c) 

In  the  general  case,  the  modulation  functions  are  more  complicated  and  involve  more  than  just  one  modulating 
frequency.  The  following  example  is  typical  of  the  modulations  we  encounter  and  is  directly  applicable  for  many 
channeled  polarimeters, 

M 

Im(x)  =  sini^mX),  (5.2) 

m= 1 

where  denotes  that  the  function  could  either  be  a  cosine  or  a  sine.  Since  we  have  M  sinusoids  multiplied 
together,  we  can  create  a  2M  x  2M  matrix  that  will  describe  all  the  possible  combinations  of  either  ±£m  of  the 
^-function,  as  well  as  distinguish  between  a  cosine  and  a  sine.  Each  sub-function  will  have  a  phasor  that,  when 
multiplied  together,  will  yield  the  net  phase  of  the  particular  channel- weight.  Employing  matrix  notation,  we 
can  create  this  “look-up-table”  by  means  of  an  outer  product  of  two  matrices: 


Rm=[Ii  fa  f m  ]  ,  (5.3a) 

£m=[pi  02  •••  O  m],  (5.3b) 


where  /m  *.  is  0  for  cosine  and  1  for  sine,  while  om  £  is  —1  for  — and  +1  for  FM  and  Ojk  are  both  2M  x  M 
in  size.  The  frequency  phase  matrix  (FPM)  is  then, 


Em  =  gar  exP 
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EmOm) 


2  v-  — 


(5.4) 


Cases  of  M  =  1, ...  ,4  can  be  seen  in  Figure  [7T7T| 

For  some  polarimeters  the  induced  modulation  may  be  more  complicated  than  the  one  prescribed  by  Equa¬ 
tion  [572]  and  the  one  directly  solved  by  creating  a  single  FPM.  In  order  to  treat  those  modulation  schemes,  we 
can  allow  for  addition  of  modulating  functions,  effectively  treating  individual  FPMs  as  a  basis  set.  The  total 
FPM  can  be  decomposed  into  L  sub-FPMs, 


E  total  —  Em  q  +  Em1  +  •  ■  ■  +  E  M  L_11 


(5.5) 


which  can  be  calculated  separately  and  simply  added  together.  The  involved  coefficients  will  be  more  complicated, 
but  they  should  present  no  additional  challenge  within  the  prescribed  methods. 
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Figure  5.1:  First  four  frequency  phase  matrices.  The  circles  represent  the  polar  form  of  the  coefficients.  Each  element 
has  an  implied  weight  of  2_JU,  which  is  omitted  to  avoid  confusion.  Also  note  that  if  ©  =  ©.  then  certain  (5-functions  will 
combine  at  5(£  ±  •■•+&  —  £7  ±  •■• )  and  (5(£  ±  •  •  •  —  £i  +  £j  ±  •  •  • ).  The  magnitudes  of  those  impulses  will  change.  The 
side  brackets  denoted  with  2/4/8/16  can  be  used  as  crop  guidelines  for  obtaining  FPMs  for  M  =  1/2/3/4.  The  circles 
correspond  in  the  following  way:  G  =  +1,  •  =  +j,  G  =  —1,  G  =  —j. 


We  have  several  ways  by  which  to  combine  multiple  dimension  modulations  into  a  total  structure  of  channels. 
In  the  case  that  we  have  determined  each  dimension’s  structure,  we  can  combine  them  using  a  Kronecker  product, 
namely, 


q{T}  ©  qM  (g>  q{4}  <g>  q{l)}.  (5.6) 

On  the  other  hand,  if  the  order  of  modulation  dimensions  is  alternating  between  elements,  we  could  use  convolu¬ 
tion  to  create  the  JV-dimensional  cloud  of  channels  that  we  would  then  unfold  into  a  vector.  As  an  example,  we 
will  consider  four  polarization  modulation  elements  that  operate  over  x/y/x/y  or  modulate  into  £/t?/£/t?.  The 
total  vector  is  then 

vec  (sUSeJ  *™  SlW  3(«e3}  q{„e4})  >  (5-7) 

where  *n  redundantly  implies  that  the  vectors  are  differently  oriented  or,  more  generally,  can  be  described  as 
degenerate  iV-dimensional  structures.  In  this  example,  q{jei}  and  q{^e3}  are  row  vectors,  while  q{,/e2}  and 
q{,)c4}  are  column  vectors.  The  result  of  the  convolution  operation  is  a  matrix  and  needs  to  be  unfolded  using 
vec  operation  that  we  defined  above.  The  choice  of  row/column  over  column /row  addressing  is  arbitrary  at  first, 
but  once  chosen  must  be  maintained  consistently  throughout  the  analysis. 

Another  way  of  generating  the  sought-after  total  vector  is  by  recognizing  that  the  Mueller  element  modulation 
patterns  can  be  alternately  viewed  as  either  a  test  dyad  or  a  projection  target.  By  treating  it  as  a  dyad, 
D  =  AGt,  we  can  look  at  its  Fourier  transform,  ©"{D}  =  T  {A}  *  T  {G}1,  with  *  now  being  a  matrix 
convolution  (same  as  multiplication,  but  every  product  is  replaced  with  convolution  between  the  same  elements 
and  added  as  before).  That  allows  us  to  combine  PSA  and  PSG  modulations  as 

qmiJ  =  vec(q9i  *  qQj)  . 


(5.8) 


Filially,  using  any  of  the  prescribed  methods  to  construct  the  Mueller  element  modulation  vectors,  we  need 
to  combine  all  of  them  into  the  corresponding  Q  matrix 
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that  maps  an  input  Mueller  vector  into  a  channel  vector, 


(5.9) 


(5.10) 


where  C  describes  the  channels  contents.  However,  we  want  to  perform  the  opposite  operation,  since  we  measure 
channels  directly.  To  do  that,  we  can  simply  take  the  pseudo  inverse  of  Q,  like  we  do  for  the  DRM  in  regular 
polarimeters, 


£+  =  (£tg)-1£t  =  (gt9)\£t.  (5.11) 

By  correctly  arranging  Fourier  transform  operations  around  the  multiplication,  we  can  use  that  reverse  mapping 
to  get  back  to  the  Mueller  elements’  information, 

M'  =  J7-1  |  Q+^r{C}|  .  (5.12) 

An  important  piece  of  insight  can  be  obtained  if  we  recognize  that  Q  is  not  that  much  different  from  D however, 
if  before  we  relied  on  constructing  multiple  dyads  against  which  to  test  the  Mueller  object,  we  can  now  have 
a  very  limited  number  of  dyads.  This  is  because  the  particular  modulation  choices  create  a  multi-dimensional 
“pointer”  that  can  be  unfolded  to  the  full  Q  representation. 

5.2  Snapshot  Channels 

A  snapshot  measurement  implies  that  there  is  no  temporal  modulation,  which  consequently  means  that  no 
additional  (other  than  the  exposure  time)  temporal  band-limit  constraints  are  placed  on  the  captured  scene.  If 
we  consider  that  each  dimension  on  the  2D  detector  can  carry  spatial,  spectral  or  no  modulation,  we  can  create 
a  verbose  set  of  nine  snapshot  channeled  systems  that  can  be  seen  in  Figure  |5.2[ 


(a)  None-None  (b)  None-Spec  (c)  None-Spat 


(d)  Spec-None  (e)  Spec-Spec  (f)  Spec-Spat 


(g)  Spat-None  (h)  Spat-Spec  (i)  Spat-Spat 


Figure  5.2:  Snapshot  systems.  Case  (a)  provides  no  modulation.  Case  (e)  is  not  straightforward  to  implement  physically. 
Case  pairs  (b)/(d),  (c)/(g)  and  (f)/(h)  are  essentially  equivalent. 


5.3  Multiple- Snapshot  Channels 

Further  developing  the  consideration  of  physical  realizability,  we  can  take  several  snapshot  measurements.  This 
gives  an  easy  access  to  a  third  modulation  dimension  —  time.  Provided  that  the  temporal  modulation  is  captured 
in  even  time  steps,  we  can  use  the  Fourier  transformation  to  create  our  channels.  However,  since  in  most  cases 
this  will  create  more  channels  than  the  original  data,  with  all  channels  being  a  single  pixel,  using  the  Fourier 
coefficients  does  not  present  any  advantage.  Instead,  the  captured  snapshots  themselves  can  be  used  as  “direct 
channels”,  or  simply,  projection  targets  like  in  the  W'  formalism,  namely, 


Q  total 


(5.13) 


This  removes  the  need  to  have  evenly  spaced  samples,  yet  maintains  the  compressed  nature  of  Q.  Note  that  even 
though  we  selected  direct  channels  only  for  temporal  modulation,  it  is  possible  to  treat  other  domain  modulations 
similarly. 


5.4  Channeled  Reconstruction 

We  use  an  SVD  method  to  calculate  the  pseudoinverse  for  reasons  of  its  numerical  stability.  First,  Q  is  decom¬ 
posed  as 

q  =  (5.i4) 

where  the  matrices  U  and  V  are  TV  x  TV  and  16  x  16  complex,  orthogonal  matrices,  respectively,  and  S  is  a 
TV  x  16  reduced  diagonal  matrix  containing  the  TV  singular  values  er i  >  cr2  >  . . .  >  <tjv  >  0.  The  pseudoinverse 
can  be  written  as 

=  (5.15) 

where  S+  is  the  16  x  TV  reduced  diagonal  matrix  containing  the  inverse  of  the  singular  values.  The  rank  of  the 
measurement  can  be  calculated  as  tr(Q  +  Q).  It  is  interesting  to  look  at  the  diagonal  elements  of  this  matrix 


to  see  how  the  information  from  the  N  measurements  is  distributed  in  the  estimated  Mueller  matrix.  We  define 
the  “reconstructables”  matrix 

B' =  vec(B)  =  diag(Q+ Q).  (5.16) 

For  each  Mueller  k- th  element,  \JWk  tells  us  the  fraction  of  energy  that  is  maintained  after  reconstruction.  When 
b'k  =  0,  sensor  space  and  scene  space  are  orthogonal,  TZ  _L  S.  Equivalently,  it  can  be  said  that  the  information 
lies  in  the  null  space  of  the  measurement.  P  When  b'k  =  1,  then  TZ  C  S  and  the  information  can  be  reconstructed 
to  within  noise  limitations. 

Although  B  can  be  useful,  it  is  a  summary  metric  and  does  not  tell  the  full  story  about  the  polarimeter.  If 
we  investigate  the  Q  +  Q  multiplication,  we  can  gain  more  insight,  since 


=  YqI+qIqYV 


(5.17) 

(5.18) 


Sq  S  q  will  reveal  a  sub-identity  matrix  that  will  have  first  R  (rank  of  Q)  elements  of  the  diagonal  equal  to 

unity  and  the  remaining  16  —  R  elements  equal  to  zero.  This  has  the  effect  of  cropping  V  into  two  —  the  column 
space  that  the  polarimeter  supports  and  the  null  space  that  the  polarimeter  rejects.  By  keeping  only  the  column 
space  vectors,  we  can  truncate  V,  namely  into  V',  that  can  be  used  to  rewrite  Equation 


5.16 


B '  =  diag  (Yq  Yq ) . 


(5.19) 


A  full  system  will  have  a  reconstructables  vector  of  all  ones,  while  partial  systems  could  design  b'k  based  on  prior 
knowledge  of  the  scene  in  question.0'  This  concept  is  developed  more  fully  for  general  pMMPs  elsewhere.18 


5.5  Examples 

By  simplifying  the  generation  of  a  polarimeter  design,  we  can  calculate  the  idealized  SNR  from  Q  directly, 
without  performing  full  simulations  of  the  system.  Simply  changing  the  way  the  problem  is  written  allows  us  to 
introduce  more  optimization  parameters  that  can  help  us  find  an  optimal  polarimeter.  In  this  section  we  will 
look  at  systems  and  see  how  the  introduced  concepts  could  help  increase  performance. 

Sabatke  introduced  Equally  Weighted  Variance  (EWV)  as  an  appropriate  metric  to  evaluate  Stokes  polarime- 
tersp2  and  Twietmeyer  later  adopted  a  similar  metric  for  use  with  Mueller  polarimetersP^ 


EWV  =  tr 


K  W/+ 


E 


15 
k- 0 


V^W'  +  .fc’ 


(5.20) 


with  K  W/+,  the  covariance  matrix,  W'+W,+T.  For  the  purposes  of  designing  pMMPs,  it  is  necessary  to 
generalize  the  combination  of  Mueller  element  variances  with  arbitrary  preference.  Instead  of  implying  a  uniform 
weight  for  all  elements,  we  can  introduce  a  weighing  vector,  u,  to  calculate  a  Weighted  Variance  (WV), 


WV  =  u1  diag 


K 


W'  + 


W'  +  ,fc- 


For  channeled  systems,  we  can  calculate  the  WV  similarly  as 


WV=uT  diag 


Kq  + 


v^lo  / 

=  Efc= 0  Uk/c 


Q  +  ,fc> 


(5.21) 


(5.22) 


with  Kq+  =  Q+Q+^,  the  product  which  contains  the  Mueller  element  reconstruction  variances  within  its 
diagonal. 


5.5.1  Spectrally  Channeled  Polarimeter 


Consider  the  spectral-none  channeled  polarimeter  proposed  by  HagenJ^  an  example  of  a  system  in  Figure  5. 2d 
The  measured  intensity  can  be  written  as 


3  3 

i^)  =  E  E  0)ma  0)f*  0)’ 

i—0  j— 0 


(5.23) 


where  •&  is  used  to  denote  a  set  of  domains  where  the  information  is  modulated.  For  this  system,  there  is  only 
modulation  in  wavenumber,  meaning  i9  =  {a}  and 


The  argument 


1 

COs(c4Cr) 

sin(c3tr)  singer)  ’ 
cos(c3cr)  sin(c4 a) 

1 

cos(ci  a) 

sin(cicr)  sin(c2<T)  ’ 
sin(ci a)  cosher) 


CiC  =  2irTicr  =  2'Kd0di\0Ba ,  r,  =  d0diX0B 


(5.24a) 


(5.24b) 


(5.25) 


contains  the  global  thickness  factor,  d0,  individual  retarder  thickness  factor,  d,,  center  wavelength,  A0,  and 
birefringence,  B.  The  vector  d  contains  all  the  modulation  information  and  will  be  our  optimization  target. 
Hagen  chooses  d  =  (  1  2  5  10  ),  with  the  resulting  channels  and  the  proposed  reconstruction  scheme. 


From  the  proposed  reconstruction,  we  see  that  some  measurements  are  ignored  for  the  sake  of  algebraic 
simplicity  —  only  channels  Co  -  Cio  are  referenced,  with  real  and  imaginary  operators  constituting  the  use  of 
conjugates.  Thus,  instead  of  using  all  37  channels,  only  21  are  used.  An  alternative  method  would  be  to  recognize 
the  modulation  induced  by  retarders,  construct  an  appropriate  FPM,  look  up  the  coefficients  and  construct  Q 
by  placing  them  at  the  contributing  frequencies.  The  resultant  Q  and  its  inverse  can  be  seen  in  Figure 


5.3 


Cn 

Channel  Content  x(64/S,in  o) 

0 

16m0o 

±i 

8m0i  ±  4m02  ±  4im03 

±2 

-m22  ±  im2 3  T  im3 2  -  m33 

±3 

-4?n02  =F  4 im03 

±4 

2to2i  ±  m2 2  =F  im23  ±  2im3i  ±  im32  ±  m33 

±5 

4m2o  ±  4  jm3o 

±6 

2m2i  ±  m22  ±  lm23  ±  2im3i  ±  im32  -  m33 

±7 

— 2m\2  ±  2im\3 

±8 

—m22  =F  im23  T  im32  +  m33 

±9 

4mn  ±  2mi2  ±  2 im\3 

±10 

O 

00 

±11 

4mn  ±  2mi2  ±  2i?ni3 

±12 

m2 2  =F  *77123  T  *m32  ~  m.33 

±13 

— 2toi2  =F  2fmi3 

±14 

— 2m2i  -  m2 2  ±  im23  ±  2im31  ±  im32  ±  m33 

±15 

-4m2o  ±  4 im30 

±16 

— 2m2\  —  m2 2  =F  im2 3  ±  2im3i  ±  im32  —  to33 

±18 

m22  ±  im2 3  =F  im3 2  +  m33 

Table  1:  Hagen  polarimeter’s  channels. 


mjjjcr) 


m00(a)  = 

4c0 

TOoi(cr)  = 

8(ci  +  c3) 

m02(a)  = 

-165R[c3] 

m03(a)  = 

165R[ci] 

toio(ct)  = 

8cio 

mn(cr)  = 

16(c7  +  c9) 

m-i  2  (cr)  = 

-325i[c7] 

mi3(cr)  = 

323[c7] 

m.20(a)  = 

165R[c5] 

m2i(cr)  = 

325R[c2  +  c4] 

m22{a)  = 

-32  K[c2  +  c8] 

m23  (ct)  = 

323  [c2  -  c8] 

m30  (ct)  = 

16SR[c5] 

to3i(ct)  = 

163[c2  +  c4  +  c6  +  c8] 

m32(a)  = 

-323[c2  +  c8] 

m33(a)  = 

325R[c8  -  c2] 

Table  2:  Hagen’s  proposed  reconstruction. 
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Figure  5.3:  Hagen’s  spectropolarimeter,  d  =  (  1  2  5  10  ).  The  matrix  containing  21  cropped  channels  can  be  seen 
between  the  two  horizontal  lines  and  has  an  EWV  of  355;  including  the  other  16  channels  lowers  the  EWV  to  187.  Note 
that  these  extra  channels  must  be  measured  in  order  to  prevent  aliasing.  The  distinction  is  whether  the  data  contained 
within  these  channels  is  used  in  reconstruction,  after  the  Fourier  transform  of  the  measured  intensity  was  determined. 
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Figure  5.4:  By  swapping  the  first  two  elements,  i.e.,  d  =  (  2  1  5 


10  ),  we  can  lower  the  EWV  down  to  130|. 


The  system  in  Figure  |5.4|  was  found  by  trial-and-error  to  see  if  any  other  arrangement  of  exactly  the  same 
elements  will  produce  better  results.  The  system  in  Figure pT4| was  found  by  optimizing  EWV  with  the  preference 
towards  fewer  channels.  Changing  to  d  =  (  2  1  4  11  )  has  the  effect  of  “orthogonalizing”  the  channels  in  a  way 
that  the  PSA-channels  are  available  independently  from  PSG-channels,  a  characteristic  that  we  have  empirically 
observed  to  be  indicative  of  optimality.  Considering  all  the  channels  in  Figure  5.3  lowers  EWV  by  47.3%,  while 
the  systematic  approach  to  measurement  selection  brings  another  35.3%  reduction  to  EWV.  In  total,  EWV  was 


reduced  to  34.1%  of  its  original  value,  suggesting  that  the  polarimeter’s  SNR  is  almost  three  times  higher.  Since 
the  number  of  channels  did  not  increase,  no  spectral  resolution  was  lost. 


5.5.2  Spatially  Channeled  Polarimeter 

Next  consider  the  spatial-spatial  channeled  polarimeter  described  by  Kudenovp  an  example  of  a  system  in 
Figure  |5.2i|  The  modulation  is  achieved  via  polarization  gratings  that  separate  the  different  Mueller  matrix 
elements  onto  patterns  of  frequencies  that  are  determined  by  the  spacing  of  the  elements.  The  intensity  can  be 


similarly  represented  as  in  Equation  5.23  with  i?  =  {x,  y}.  Kudenov  used  the  following  modulations: 


f  A(x,y) 


lc{x,y ) 


1 

cos(47ra:) 

sin(47ra)  cos(47tj/)  ’ 
sin(47ra;)  sin(47ry) 

1 

cos(27 iy) 

sin(27ry)  cos(27ra:) 
sin(27rj/)  sin(27ra:) 


(5.26a) 


(5.26b) 
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(c)  &x  -  Ciy  - 

(d)  ax  = 

2gx  =  2  gy  —  2, 

2gx  —  2  gy 

x/y/y/x,  33/49 

x/y/x/y, 

channels  used, 

channels 

EWV  =  209 

EWV  =  15: 

ay  =  (e)  a i  =  2a2  = 

=  2,  gi  =  2g2  =  2, 

35/49  x/x/y/y,  49/49 

used,  channels  used, 

EWV  =  121 


Figure  5.5:  Top  row  shows  the  Ov  plane  of  channels  of  each  configuration.  Bottom  row  shows  Kq  +  . 


Figure  [575] shows  a  comparison  between  three  systems  with  the  only  difference  being  the  order  of  modulation. 
The  merit  of  introducing  Q  is  clear;  better  performance  is  achieved  virtually  for  free,  using  the  same  polarization 
elements  in  a  different  order.  Systems  in  Figure  |5.5|  were  found  by  optimizing  using  genetic  algorithm,  while 
continuously  relaxing  the  design  restrictions. 
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Figure  5.6:  Frequency  grid  of  the  Mueller  modulation.  £  and  r/  are  the  and  y-axes,  respectively. 
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Figure  5.7:  Optimal  Spatial-Spatial  Polarimeter. 


Although  a  symmetrical  (x/y/y/x)  modulation  design  may  seem  intuitively  logical,  it  is  possible  to  improve 


the  design  as  evidenced  by  the  polarimeter  in  Figure  |5.5e|  which  is  shown  in  greater  detail  in  Figure  |5.6|  and 

an  asymmetrical  order  of  modulations  improves  EWV  by  27.8%.  Then,  in 


First,  in  Figure  5.5d 


Figure 

Figure  5.5e|  it  is  improved  by  another  19.9%  by  splitting  the  modulation  into  one  dimensional  structures  for 
the  PSG  and  the  PSA.  In  total,  EWV  was  reduced  to  57.9%  of  its  original  value.  Although  not  as  large  of  an 
improvement  as  in  the  previous  example,  it  is,  nonetheless,  significant. 


The  reason  for  the  EWV  improvements  lies  in  how  the  channels  interfere.  From  the  comparison  of  Figure  [53] 
we  can  note  that  the  better  systems  “focus”  the  reconstruction  into  the  diagonal  of  K  q  .  Upon  further  investiga¬ 
tion,  it  is  possible  to  make  a  general  statement  that  we  want  Q  to  be  a  matrix  that  is  unitary  through  a  scalar. 
When  this  condition  holds,  the  channel  structures  created  for  each  of  the  Mueller  elements  forms  an  orthogonal 
basis  set,  which  will  force  all  the  variances  to  lie  on  the  diagonal,  thereby  ensuring  the  minimum  achievable  EWV 
for  the  number  of  modulations  introduced.  In  fact,  we  can  calculate  the  minimum  EWV  simply  as 


EWVmin  =  ^nG®  n^, 


(5.27) 


where  n  is  a  vector  that  contains  the  number  of  channels  that  the  PSG’s  and  the  PSA’s  Stokes  parameters  are 
split  into.  If  we  are  to  assume  a  spherical  modulation,  i.e., 


£  A 


f  G 


1 

cos(27ra£i) 

sin(27ra;^i)  cos(27ra"^2)  ’ 

sin(27ra;£i)  sin(27ra:^2) 

1 

cos(27rj/r7i) 

sin(27r?/?7i)  cos(27ry?72)  ’ 

sin(27ry?7i)  sin(27r(/772) 


(5.28a) 


(5.28b) 


then  nj  =  nn  =  1  2  4  4 


and  the  minimum  EWV  is  121,  like  in  Figure 

iT 


5.5e 


A  better  EWV  is  math- 
in  which  case  EWVmin  =  81.  However,  there 


ematically  possible  if  for  example,  =  nc  =  [  1  2  2  4 

appears  to  exist  no  such  physically  realizable  modulation  scheme  as  it  would  require  analyzing  and  generating 
vectors  to  reach  a  DoP  of  a/2- 


5.5.3  Multiple  Snapshot  Polarimeter 

As  mentioned  before,  our  matrix  Q  is  nothing  more  than  a  dyad  that  contains  modulations  within  each  Mueller 
channel  projection  target.  Thus,  once  we  determine  the  system’s  spectral  and/or  spatial  modulations,  we  can 
rotate  the  channel  structures  via  a  unitary  transformation  that  will  remix  the  channel  structures  in  PSG  and 
PSA 


Qe  = 


UaJ-{A} 


(5.29) 


Applying  Ug  and  will  keep  the  channels’  relative  orientation,  thereby  maintaining  the  EWV.  Looking 
at  the  Mueller  matrix  of  a  linear  retarder^!  We  can  show  that  |  det(LR(<5,  0))|  =  1,  which  means  a  unitary 
transformation.  Thus,  we  can  construct  a  multiple  snapshot  system  simply  by  enclosing  the  system  between  two 


retarders  that  have  different  orientations  for  each  measurement.  Using  Equation  5.13  gives  us  the  total  Q. 
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(a)  One  snapshot  —  49  channels 
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(b)  Two  snapshots  —  98  channels 
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(c)  Three  snapshots  —  147  channels 
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(d)  Four  snapshots  —  196  channels 


Table  3:  Optimization  results  for  different  number  of  temporal  snapshots 


We  ran  optimizations  of  64  differently  configured  spatial-spatial  channeled  polarimeters.  The  distinction  is 
in  the  frequencies  used  for  Oi,  a 2,  g\,  g2  and  the  corresponding  dimension  into  which  the  data  was  mapped. 
Retardances  and  orientations  of  the  retarders  were  the  other  variables.  A  genetic  algorithm  was  used  to  find  the 
lowest  EWV. 

From  these  results  we  gather  that  as  the  number  of  temporal  measurements  grows,  the  importance  of  the 
spatial  frequencies  and  order  of  modulations  diminishes.  This  bodes  well  if  we  are  to  understand  this  phenomenon 
as  a  continuously  growing  temporal  bandwidth  constraint  allowing  us  to  simplify  the  spatial  multiplexing. 

5.6  Conclusion 

Introducing  Q  and  methods  for  generating  it  automatically  allows  us  to  describe  a  wide  range  of  similar  systems 
with  a  handful  of  parameters  and  removes  the  need  to  handle  reconstruction  by  hand.  Furthermore,  analysis  of 
Q  reveals  certain  design  metrics  immediately  instead  of  having  to  run  an  elaborate  simulation.  The  end  result  is 
that  a  more  optimal  system  can  be  found  often  without  requiring  the  use  of  any  extra  elements,  while  injecting 
the  optimization  procedure  before  element  selection  will  allow  us  to  choose  better  elements. 


6.  PARTIAL  POLARIMETERS 


Numerous  authors  have  studied  the  structure  of  the  Mueller  matrix,  and  much  is  known  about  how  the  various 
Mueller  matrix  elements  relate  to  the  physical  properties  of  diattenuation,  retardance,  and  depolarization.31  °9 
It  should  be  clear  that  not  all  4  x  4  real  matrices  are  physically  realizable.  A  physical  Mueller  matrix  must  map 
real  sets  of  Stokes  parameters  into  real  sets  of  Stokes  parameters,  but  there  are  other  conditions  that  must  also 
be  met  as  recently  discussed  by  Gil.® 

Much  of  the  literature  on  Mueller  matrices  is  concerned  with  methods  to  decompose  the  Mueller  matrix  in 
order  to  understand  its  structure  and  relate  it  to  scattering  properties.  In  the  class  of  series  decompositions, 
the  Mueller  matrix  is  broken  up  into  discrete  diattenuating,  retarding,  and  depolarizing  layers,  and  the  result 
is  a  product  of  Mueller  matrices  that  describe  the  effects  of  the  whole.  Lu  and  ChipmaiJ®  developed  a  series 
decomposition  that  writes  the  Mueller  matrix  as  a  non-unique  cascade  of  pure  diattenuation,  retardance,  and 
depolarization  Mueller  matrices.  Ossikovski  and  colleagues  developed  a  different  decomposition  that  eliminated 
the  order-dependence  of  the  Lu-Chipman  decomposition  by  creating  a  decomposition  that  is  symmetric  through 
the  Minkowski  metric  tensor  G  =  diag(l,  —  1,  —  1,  —  1).®  It’s  clear  that  while  one  can  use  either  of  these 
decompositions  (or  any  other),  they  may  not  actually  represent  the  physics  of  any  particular  process. 

The  limit  of  series  decompositions  is  the  class  of  differential  decompositions.^2  These  split  the  Mueller  matrix 
into  differential  slices  in  an  attempt  to  identify  its  fundamental  characteristics.  Noble  and  Chipman31  63  use  the 
method  of  matrix  roots  to  uncover  a  fundamental  differential  Mueller  matrix  that  can  be  written  in  terms  of  15 
Mueller  matrix  generators  -  three  for  retardance,  three  for  diattenuation,  and  nine  for  depolarization.  Ossikovski 
developed  a  logarithmic  decomposition  of  the  Mueller  matri:x®  that  operates  using  a  different  formalism,  but 
produces  an  equivalent  outcome  to  that  of  the  matrix  roots  decomposition.® 

A  third  class  of  decomposition  is  the  class  of  additive  decompositions  that  consider  the  Mueller  matrix  as 
an  ensemble  average  of  parallel  scattering  processes  that  are  added  incoherently.  Gil  provides  a  recent  review 
that  covers  the  general  cases  of  the  trivial,  spectral,  and  arbitrary  decompositions.00  The  most  famous  parallel 
decomposition  is  that  of  Cloude,  who  demonstrated  that  an  arbitrary  Mueller  matrix  could  be  written  as  a 
superposition  of  not  more  than  four  pure  Mueller-Jones  matrices^65  Ossikovski  has  demonstrated  rigorously  that 
in  the  limit  of  weakly  depolarizing  Mueller  matrices,  all  decompositions  return  identical  polarization  properties 
to  first  order.®  However,  for  more  general  depolarizing  matrices,  the  various  methods  return  different  results 
for  the  “fundamental”  properties  or  retardance  and  diattenuation  of  a  Mueller  matrix  under  test. 

All  of  these  classes  of  decompositions  are  important  for  understanding  the  fundamental  properties  of  the 
Mueller  matrix.  However,  measurement  of  the  Mueller  matrix  requires  us  to  consider  a  different  basis  set 
altogether.  A  Mueller  matrix  polarimeter  operates  by  using  a  polarization  state  generator  (PSG)  to  illuminate 
the  sample  with  a  controlled  state  of  polarization.  The  polarimeter  then  measures  the  irradiance  passed  through  a 
polarization  state  analyzer  (PSA)  set  to  a  second  polarization  state.  Through  a  suitably  diverse  set  of  illumination 
and  analysis  states,  the  elements  of  the  Mueller  matrix  can  be  determined.®  Much  as  is  the  case  in  Stokes 
polarimetry,®  the  measurement  corresponding  to  each  pair  of  PSG/PSA  states  can  be  thought  of  as  a  projection 
onto  a  basis  vector,  and  then  the  unknown  Mueller  matrix  can  be  estimated  through  a  least-squares  inversion 
process  that  produces  an  additive  decomposition.  Once  the  problem  is  cast  in  this  manner,  the  design  of  a 
measurement  system  then  becomes  an  optimization  problem  where  a  particular  measurement  basis  is  chosen 
in  order  to  highlight  specific  aspects  of  the  Mueller  matrix.  At  least  16  measurements  are  needed  in  order  to 
reconstruct  the  full  Mueller  matrix  in  general]*-8  while  the  choice  of  specific  illumination  states  can  help  balance 
the  signal-to-noise  ratio  (SNR)  and/or  error  in  particular  Mueller  matrix  elements.'08  66  Going  one  step  further, 
we  can  even  design  a  partial  Mueller  matrix  polarimeter  (pMMP)  that  allows  certain  elements  or  combinations 
of  elements  of  the  Mueller  matrix  to  be  recovered  with  fewer  than  16  measurements  while  ignoring  other  elements 
that  might  not  be  necessary  for  a  particular  sensing  problem.®  Hoover  and  his  coworkerdUM  have  demonstrated 
that  reduced  dimensionality  subspaces  of  Mueller  matrix  space  can  be  used  to  perform  invariant  target  detection 
through  nonlinear  model  fitting.  Goudail  and  his  coworkers6'  69  have  demonstrated  that  a  single-measurement 
pMMP  is  optimal  for  maximizing  polarization  contrast  in  a  two-class  detection  problem  with  known  class  Mueller 
matrices. 

In  this  paper  we  consider  the  design  of  pMMPs  that  seek  to  measure  certain  aspects  of  the  Mueller  matrix 
that  might  be  dictated  by  a  particular  sensing  task.  The  pMMP  could  be  an  imaging  or  non-imaging  device,  but 


the  design  of  the  instrument  proceeds  from  knowledge  of  linear  combinations  of  Mueller  matrix  elements  that 
allow  a  particular  task  to  be  performed.19  20  It  is  well  known  that  it  is  not  possible  to  measure  a  single  Mueller 
matrix  element  or  a  single  arbitrary  combination  of  Mueller  matrix  elements  in  a  single  measurement  due  to 
the  restrictions  on  the  structure  of  the  Stokes  parameters  of  the  PSG  and  PSA  states.  Previous  author£5l'Pl 
have  considered  specific  optimization  strategies  designed  to  maximize  performance  on  a  particular  task.  In  this 
paper  we  approach  the  more  general  two-part  problem  of  a)  identifying  the  proper  subspace  in  which  to  make  a 
detection  decision  and  b)  designing  a  pMMP  to  get  as  close  as  possible  to  a  specified  subspace  of  Mueller  matrix 
space  through  careful  selection  of  measurement  states.  In  order  to  accomplish  this,  we  discuss  some  of  the  details 
of  the  structure  of  the  Mueller  matrix  and  how  it  interacts  with  the  PSG  and  PSA  before  developing  a  numerical 
optimization  method  that  produces  the  desired  pMMP  design. 


The  remainder  of  this  manuscript  is  organized  as  follows.  Section  |6.1|  reviews  the  mathematics  of  Mueller 
matrix  polarimetry  and  discusses  the  modifications  necessary  to  consider  pMMPs.  Section  |6.2|  considers  the 
structure  of  a  few  pMMPs  in  a  way  that  allows  us  to  understand  how  the  PSG  and  PSA  interact  with  the 
Mueller  matrix  to  build  up  a  pMMP  basis.  In  Section  |63|  we  generalize  patterns  seen  in  Section  [612] to  a  general 
class  of  4 ij  pMMP  systems,  as  well  as  develop  various  metrics  by  which  to  evaluate  the  noise  resilience  and  the 
proximity  of  a  Af-dimensional  subspace  of  Mueller  matrix  space  to  an  Admeasurement  pMMP.  In  Section  |6.4| 
we  apply  the  developed  concepts  to  an  object  discrimination  task  from  the  literature19  and  discuss  the  results. 
Section  6.5  concludes  the  paper. 


6.1  Notation  Overview 

Before  discussing  the  topics  of  pMMP  design,  we  must  first  review  the  concepts  of  polarimetric  sensing  and 
expand  them  in  a  way  so  as  to  expose  its  structure  and  thereby  give  us  tools  to  handle  the  information  mapping 
successfully. 

Mueller  Matrix  Polarimetry 

Consider  an  unknown  Mueller  matrix  M  that  modifies  the  Stokes  parameters.  The  system  PSG  generates  an 
incident  beam  with  Stokes  parameters  G  and  the  PSA  analyzes  the  scattered  light  by  taking  a  projection  of  the 
scattered  Stokes  parameters  onto  the  state  described  by  the  Stokes  parameter  set  A.  The  nth  measurement  of 
irradiance  in  the  polarimeter  is 


4  =  aImg„  =  d'"m'.  (6.1) 

In  Eq.  |6.1|  we  define  the  16  x  1  vectors: 

Dn  =  vec(  A„G^)  =  An  ®  G„,  (6-2) 


and 

M'=vec(M),  (6-3) 

where  (g)  is  the  Kronecker  (direct)  product,  and  vec(M)  creates  a  column  vector  by  reordering  the  matrix  M 
into  a  vector  in  a  row- by- row  fashion.  Eq.  |6.1|shows  that  a  single  measurement  of  the  pMMP  is  a  projection  of 
the  unknown  Mueller  matrix  onto  a  known  basis  vector  in  R16.  By  taking  a  collection  of  such  projections,  the 
unknown  matrix  -  or  portions  of  it  in  the  case  of  a  pMMP  -  can  be  determined  in  a  least-squares  sense.  The 
series  of  N  measurements  in  a  polarimeter  is 
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+  n  =  W  M'  +  n. 


In  Eq.  6.4  the  TV  x  16  measurement  matrix  W  is 


W  =  {  D(  D' 


(6.4) 


(6.5) 


and  n  is  considered  to  be  additive  noise  with  variance  w2.  The  Mueller  matrix  estimation  is  typically  accom¬ 
plished  by  using  the  pseudoinverse: 


m'  =  W  +  I  =  W  + WM'  +  w+n. 


(6.6) 


6.6 


The  performance  of  M  in  Eq. 
pseudoinverse  operates  on  the  noise. 


is  dictated  by  how  close  W  +  W  is  to  the  identity  matrix  and  how  the 


Structure  of  Partial  Mueller  Matrices 

In  the  case  of  partial  polarimeters,  N  <  16,  and  the  maximum  rank  that  the  polarimeter  can  achieve  is  N.  It  is 
easy  to  demonstrate 

tr  (^  +  ^)  =  rank  (^)  =  R.  (6.7) 

In  this  paper  we  will  use  the  singular  value  decomposition  (SVD)'1  to  compute  the  pseudoinverse.  The  SVD  of 
W  yields 

W  =  US  VT.  (6.8) 


The  matrices  U  and  V  are  Rx  R  and  16  x  16  real,  orthogonal  matrices,  respectively,  and  S  is  a  R  x  16  reduced 
diagonal  matrix  containing  the  R  singular  values  cti  >  02  >  •  •  ■  >  &r  >  0.  The  columns  of  U  span  the  range 
of  the  pseudoinverse,  and  the  columns  of  V  span  Mueller  matrix  space.  The  first  R  columns  correspond  to  the 
non-zero  singular  values  and  span  the  portion  of  Mueller  matrix  space  that  the  pMMP  can  reconstruct.  The 
pseudoinverse  can  be  written  as 

yy  +  =  Y|(  +  IZT,  (6.9) 

where  S+  is  the  16  x  R  reduced  diagonal  matrix  containing  the  inverse  of  the  singular  values.  The  SVD 
pseudoinverse  creates  a  “maximally  orthogonal”  inverse. 

Examining  the  diagonal  elements  of  W  +  W  matrix  tells  how  the  information  from  the  N  measurements 
contributes  to  the  rank  and  how  that  information  is  distributed  in  the  estimated  Mueller  matrix.  Define  the 
reconstruc.tables  matrix  as: 

B'  =  vec(]§)  =  diag(^+^)  .  (6.10) 


We  will  consider  examples  of  this  matrix  in  subsequent  subsections,  but  at  this  point  we  can  say  that  B  relates 
the  percentage  of  each  Mueller  matrix  element  that  is  reconstructed  in  the  pMMP.  In  the  limit  of  N  =  16,  the 
pMMP  becomes  a  full  polarimeter,  W  +  W  =  I16xl6,  and  B  is  a  4  x  4  matrix  of  all  ones;  all  elements  of  the 
Mueller  matrix  can  be  reconstructed. 


To  understand  the  function  of  the  pMMP,  consider  the  multiplication  of  the  matrix  and  its  pseudoinverse 


yg+W  =  YJQ+II+n|iY+  =  YY+YY+-  (6.11) 

The  matrix  S  +  S  is  diagonal  with  the  first  R  elements  equaling  unity  and  the  last  16  —  R  equaling  zero.  This 
permits  the  claim  made  in  Equation  |6.7|  Thus, 

yg+W  =  Y'Y't,  (6.12) 

where  V'  is  the  16  x  R  matrix  composed  of  the  first  R  columns  of  V.  Another  way  of  interpreting  the  SVD  of 
W  is  that  V'  forms  an  orthogonal  basis  that  spans  the  subspace  of  Mueller  matrix  space  that  forms  the  domain 
of  the  particular  pMMP  represented  by  W.  Likewise,  the  columns  of  V  discarded  by  the  SVD  (corresponding 
to  singular  values  of  zero)  span  the  null  space  of  the  particular  pMMP.  However,  as  shown  later,  knowledge  of 
the  domain  alone  is  not  sufficient  to  predict  performance,  as  the  conditioning  of  the  matrix  W  is  important  in 
the  presence  of  noise  and  error. 


6.2  Examples  of  Partial  Mueller  Matrix  Polarimeters 

This  analysis  is  restricted  to  pMMPs  that  use  fully  polarized  PSG  and  PSA  states.  Goudail  and  Tycf69  demon¬ 
strated  that  partially  polarized  PSG  and  PSA  states  never  improve  contrast.  Below  we  will  consider  cases  where 
one  or  more  PSG  or  PSA  state  is  unpolarized,  allowing  reconstruction  of  particular  elements  of  the  Mueller 
matrix  with  fewer  measurements  than  would  be  necessary  if  all  PSG  and  PSA  states  were  fully  polarized. 

Canonical  4~ Measurement  pMMP 

Consider  the  simple  N  =  4  measurement  pMMP  that  measures  the  co-polarized  and  cross-polarized  return  for 
both  vertically  and  horizontally  polarized  illumination.  For  compactness,  we  introduce  the  following  notation 
for  the  analyzer  and  generator  matrices 


A  =>  \  [  t  t  ]  ,  (6.13a) 

£  =>>  [  -  t  ->■  t  ]  ,  (6.13b) 


where  -*■  =  [  1 

[1-10  0 
definitions  of  the  analyzing  vector  is  needed  for  rigor 


0  0  ]  is  the  set  of  Stokes  parameters  for  ideally  horizontally  polarized  light  and  t  = 

is  the  set  of  Stokes  parameters  for  ideally  vertically  polarized  light.  The  presence  of  |  in  the 

the  polarization  sensing  systems  in  consideration  dismiss 
half  of  the  light  if  the  input  is  unpolarized.  The  set  of  four  PSG/PSA  pairs  in  Eqs.  6.13a  and  6.13b  results  in 
the  instrument  matrix 


1  10  0  1 
1  100-1 
1-100  1 
1-10  0-1 


1 

-1 

-1 

1 


0  0 
0  0 
0  0 
0  0 


0  0 
0  0 
0  0 
0  0 


0  0 
0  0 
0  0 
0  0 


0  0 
0  0 
0  0 
0  0 


0  0 
0  0 
0  0 
0  0 


(6.14) 


and  the  reconst ructables  matrix 


110  0 
110  0 
0  0  0  0 
0  0  0  0 


(6.15) 


This  is  the  well  known  result  that  four  measurements  are  needed  to  reconstruct  four  Mueller  matrix  elements, 
and  that  these  four  elements  must  come  in  a  “block”  pattern  within  the  Mueller  matrix A  similar  polarimeter 
could  be  obtained  with  a  4-measurment  combination  of  any  two  of  the  six  canonical  states  -»■  ,  t  ,  X  ,  x  ,  Q 
,  Q  ,  where  X  and  X  represent  45°  and  —45°,  and  Q  and  Q  represent  right-  and  left  circular  polarization, 
respectively. 


While  this  well-known  result  tells  how  to  design  a  pMMP  to  reconstruct  one  of  these  groupings  of  four 
elements,  it  is  not  obvious  how  to  add  additional  measurements  to  reconstruct  additional  elements  or  how  to 
design  a  pMMP  to  reconstruct  linear  combinations  of  elements  rather  than  isolated  elements  within  the  pMMP. 

Diagonal  Depolarization  Elements 

Depolarization  is  a  rich  physical  process  that  contains  significant  information  about  the  random  scattering 
properties!12  Noble  and  ChipmarEillkjll  recently  described  the  nine  degrees  of  freedom  for  depolarization.  Three 
of  these  correspond  to  randomness  in  the  diattenuation  properties  of  the  Mueller  matrix,  three  to  randomness 
in  the  retardance  properties  of  the  Mueller  matrix,  and  three  to  “diagonal  depolarization,”  which  is  related  to 
randomness  in  geometric  transformations  as  would  happen  in  multiple  scattering  or  rough  surface  scattering 
processes.  Often,  the  diagonal  depolarization  elements  are  important  for  discrimination  in  both  optical  and 
radar  tasks. ' 3 


Subsection  6.2  shows  that  each  canonical  four- measurement  polarimeter  provides  one  diagonal  element  (in 
addition  to  moo,  which  is  involved  in  all  pMMPs).  One  obvious  way  to  reconstruct  the  diagonal  elements  then 
would  be  to  use  a  12-measurement  pMMP  defined  by  the  analyzer  and  generator  matrices 


-  t  t  x  x  x  x  Q  Q  Q  Q], 

G  =>  [  -*■  t  *  t  X  x  x  X  Q  Q  Q  Q  ]  , 


(6.16a) 

(6.16b) 


which  produces  the  following  reconstruc.tables  matrix 


B  = 


1111 
110  0 
10  10 
10  0  1 


(6.17) 


In  addition  to  the  desired  diagonal  elements,  the  diattenuation  and  polarizance  vectors®  are  also  measured.  This 
12-measurement  pMMP  only  reconstructs  10  Mueller  matrix  elements,  since  each  of  the  three  canonical  pMMPs 
redundantly  reconstructs  moo- 

We  can  address  this  redundancy  by  eliminating  one  of  the  cross-polarized  measurements  in  two  of  the  canonical 
pMMPs  so  that 


A  =  5  [  ^  ^  ^  X  ^  Q  2  2  ]  ,  (6.18a) 

£  =  [  -  t  -  t  X  Q  Q  Q  ]  ,  (6.18b) 

which  produces  the  same  reconstructables  matrix  as  Eq.  |6.17[  In  this  case  the  elimination  of  redundancy  allows 
10  Mueller  matrix  elements  to  be  reconstructed  from  10  measurements. 

The  matrix  of  Eq.  |6.17|still  unnecessarily  reconstructs  the  first  column  and  the  first  row  of  the  Mueller  matrix. 
The  number  of  measurements  can  be  lowered  by  further  reducing  two  of  the  three  canonical  4-measurement 
pMMPs  to  two-measurement  pMMPs  that  make  co-polarized  measurements  only,  e.  g. 

A  =>  §  [  "*■  t  t  Q  2  ] , 

£  =>  [  -  t  -  t  *  x  Q  2  ] , 

which  produces  the  reconstructables  matrix 

T  l  l  I  I  I 

2  2 

_  _  110  0 

=  -  |  0  1  0  ■ 

_  I  0  0  1  _ 


(6.19a) 

(6.19b) 


(6.20) 


Examination  of  V  can  help  to  determine  how  the  elements  of  B  correspond  to  reconstructed  Mueller  matrix 
channels  as  discussed  in  section  |63|  This  pMMP  can  reconstruct  the  diagonal  elements  moo ;  mi  l ;  m22 ,  TO33  as 
well  as  the  elements  mio  and  moi-  In  addition  to  these  individual  elements,  the  polarimeter  can  also  reconstruct 
the  linear  combination  channels  (m2o  +  1^02) /V%  and  (r?t3o  +  rn os)/v/2.  Note  that  existence  of  reconstruction 
channels  does  not  guarantee  that  these  channels  will  have  acceptable  SNR.  These  items  are  discussed  in  greater 
detail  below. 

The  polarimeter  described  by  Eq.  |6.19a|  and  Eq.  |6.19b|  is  the  lowest  dimensionality  that  we  have  been  able 
to  find  that  reconstructs  all  three  of  the  diagonal  elements  with  fully  polarized  analyzer  and  generator  states. 
However,  use  of  unpolarized  measurements  adds  another  degree  of  freedom  and  provides  capacity  for  fewer 
measurements.  Consider  a  system  that  makes  six  canonical,  co-polarized  measurements  and  one  unpolarized 
measurement 


A=>|[0^  t  *  x  Q  2  ] , 

£  =>  [  O  -  t  *  n  Q  2  ]  , 

where  O  =  [  1  0  0  0  ]T.  The  reconstructables  matrix  for  this  polarimeter  is 


B  = 


x  2  2  2 
i  1  0  0 

|  0  1  0 

1  0  0  1 


(6.21a) 

(6.21b) 


(6.22) 


The  addition  of  the  one  unpolarized  measurement  allows  the  moo  term  to  be  reconstructed  directly,  obviating 
the  need  for  the  cross-polarized  measurements  indicated  in  Eq.  |6.19a|and  Eq.  |6.19b| 


6.3  Partial  Mueller  Matrix  Polarimeter  Design 

Based  on  the  understanding  gained  in  Section  |6T2|  we  can  remove  some  prior  constraints,  and  generalize  pMMP 
systems  to  have  arbitrary  fully  polarized  PSG  and  PSA  states  within  the  requirement  that  N  =  R. 

Interpreting  the  Reconstructables  Matrix 

It  is  unclear  from  B  alone,  which  Mueller  elements  are  grouped  together  into  combinatorial  channels.  When  a 
particular  element  of  B  is  unity,  then  we  can  reconstruct  the  corresponding  Mueller  matrix  element.  But  when 
it  is  other  than  unity,  the  element  must  appear  in  combination  with  other  Mueller  matrix  elements.  However,  the 
fact  that  B  is  derived  from  the  columns  of  V'  provides  insight  into  the  overall  subspace  spanned  by  the  pMMP. 
If  two  or  more  columns  of  V',  say  v„  and  vm  correspond  to  equal  singular  values  crn  =  arn .  then  they  span  a 
hyperplane  with  identical  geometrical  characteristics  in  the  context  of  W.  In  this  case,  any  set  of  orthogonal 
basis  vectors  in  that  hyperplane  can  be  used,  allowing  for  a  more  intuitive  grouping  of  Mueller  matrix  elements 
if  desired. 


All  of  the  previous  examples  featured  four-block  measurements.  Prior  work  was  already  familiar  with 
TOqo /rriiQ /niQ.j /rnl3  measurement,  but  a  similar  reconstruction  exists  for  off-Mueller  grid  measurement.  Con¬ 
sider  the  arbitrary  analyzing  and  generating  vectors: 


A±  =  [  1  ±a\  ±a2  ±a3  ]  , 
G±  =  [  1  ±g\  ±.92  ±93  ]  • 


(6.23a) 

(6.23b) 


A  four-block  polarimeter  would  go  through  the  following  four  combinations  of  measurements: 

,T 


A4 

g4 


5  [  A+  A+  A-  A- 

f  G+  G_  G+  G 


The  resulting  measurement  matrix  is 


W4  = 


(A+< 
(A+i 
(A- 1 
(A-  1 


G  +  )T 
G-)T 
G+)T 
G-)T 


(6.24a) 

(6.24b) 

(6.25) 


The  SVD  of  W4  in  Eq.  6.25  will  have  four  column-space  vectors  with  four  identical  singular  values.  If  a  particular 


representation  of  V  is  chosen,  then  there  is  only  one  U  to  go  along  with  it.  A  linear  combination  of  these  vectors 
corresponds  to  rotation  of  the  underlying  vectors.  This  operation  does  not  alter  the  space,  but  merely  rotates 
the  axes  by  which  that  space  is  described.  We  will  show  that  we  can  write  down  a  particular  decomposition  that 
is  relatively  easy  to  treat  analytically.  When  faced  with  more  complex  V  matrices  that  have  non-equal  singular 
values,  adding  and  subtracting  the  underlying  vectors  is  also  possible,  but  special  care  needs  to  be  taken. 

For  the  four-block  polarimeter,  the  column  space  can  be  described  with  the  following  set  of  vectors 


Y4  = 


10  0  0 
0  9i  0  0 

0  92  0  0 

0  93  0  0 

0  0  oi  0 

0  0  0  ai9i 

0  0  0  0192 

0  0  0  oi93 

0  0  a2  0 

0  0  0  a29i 

0  0  0  0292 
0  0  0  0293 

0  0  03  0 

0  0  0  0391 

0  0  0  a3g2 

0  0  0  a3g3 


(6.26) 


Figure  6.1:  Possible  sets  of  measurements  that  maintain  the  optimal  N  =  R. 
that  correspond  to  four  unity  singular  values  and 
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The  corresponding  reconst ructables  matrix  is 
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(6.28) 


Additional  Measurements 

In  order  to  expand  the  space  coverage,  we  need  to  add  more  measurements.  To  do  so  while  keeping  system  rank 
equal  to  the  number  of  measurements  means  that  each  additional  column  in  V'  needs  to  be  orthogonal  to  every 
pre-existing  one.  If  we  are  limited  by  the  assumption  of  fully-polarized  measurements,  then  the  new  analyzing 
vector  pair  A2,±  needs  to  be  orthogonal  to  the  pre-existing  analyzing  vector  pair  Ai  ±  in  the  Poincare  sphere 
space.  Thus,  once  the  first  pair  is  selected,  the  new  pair  is  bound  to  the  space  of  the  orthogonal  circle.  Once 
A.2  ±  is  chosen,  there  is  only  one  more  orthogonal  set  of  vectors  A3j=t-  that  can  be  added. 

It  is  important  to  make  the  connection  between  this  general  case  and  the  one  discussed  in  Subsection  4.B.  If 
we  make  four  measurements  from  each  combination  in  Ai  j-,  A,2.±  and  A3 ,±,  then  we  will  end  up  with  TV  =  12 
and  R  =  10.  This  is  because  each  block  is  capable  of  reconstructing  moo  011  its  own  and  measuring  it  three 
times  will  have  the  effect  of  averaging,  and  thereby  lowering  the  noise  in  its  reconstruction.  In  order  to  keep 
N  =  R  we  must  make  only  one  four-measurement  set  specified  by  Aj  ±  and  up  to  two  additional  fewer-than-four 
measurements  specified  by  A2,±  and  A3i±.  This  produces  16  possible  measurement  schemes.  By  denoting  the 
set  as  4 ij  and  requiring  that  4  >  i  >  j,  we  can  ignore  the  six  redundant  schemes  as  can  be  seen  in  Fig.  |6.1| 

Purely  for  purposes  of  simplifying  the  notation,  we  define  the  analyzing  and  generating  vectors  of  the  one-, 
two-  and  three-measurement  cases  as: 


As 

=►§[ 

A+ 

A-  A-]t, 

(6.29a) 

£3 

=*  [ 

G  + 

G+  G  ]T, 

(6.29b) 

=  2 

A+ 

A-  ]T, 

(6.29c) 

£2 

=>  [ 

G  + 

G-  ]T, 

(6.29d) 

Ai 

=►*[ 

A+ 

]T: 

(6.29e) 

£1 

=>■  [ 

G  + 

]T- 

(6.29f) 

It  can  be  shown  that  this  selection  considers  all  possible  combinations.  The  ±  only  denotes  the  operation  on 
the  ai/gi,  02/52,  and  03/53,  but  a  selection  of  a  different  vector  can  effectively  construct  all  other  combinations 
within  the  syntax  implied  by  ±. 

The  denoted  characteristics  of  these  measurements  are  only  correct  if  moo  is  known.  Thus,  they  can  only  be 
used  as  additional  measurements  and  their  measurement  sub-matrices  are: 


and  the  sub-reconstructables  matrices  are 


B3  = 


i2=2 


— 1  =  3 


The  total  measurement  matrix  of  a  4 ij  system  is 
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4  ij  = 

[ 

wj 

wj  W] 

dT. 

(6.30a) 

(6.30b) 

(6.30c) 

(6.31a) 

(6.31b) 

(6.31c) 


(6.32) 


The  constraints  placed  on  A  2. :  and  A.3^  mean  that  the  reconstructables  matrix  is  the  sum  of  the  sub-matrices 


^_4ij  —  Mt  +  Mi  +  Mi- 


(6.33) 


Structured  Decomposition 

As  before,  we  can  perform  SVD  on  the  matrix  to  find  the  space  coverage  and  noise  resilience  of  any  given 
polarimeter.  However,  in  the  case  of  being  limited  to  the  defined  class  of  4 ij  pMMP  systems,  we  can  introduce 
a  structured  decomposition 

4 jj  =  U s ,4j j  S g VjAi1,  (6.34) 

where  s  differentiates  this  decomposition  from  the  typical  SVD.  The  goal  of  this  decomposition  is  to  be  easily 
parsable  by  a  human  and  provide  an  intuitive  view  of  pMMP  properties.  The  following  are  the  structured 
matrices  for  any  4 ij  system: 
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=  \J\  diag(  N  <3  Q  <j  ), 

Y sMj  =  [  Y4  YJ  Yj  ] , 


(6.35) 

(6.36) 

(6.37) 


where  the  left  structured  sub-matrices  are: 
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the  effectively  rotated  singular  values  are: 


?4  —  {4, 4, 4}  , 
=  {3,3,3}  , 
7  =  {4,2}, 

«Ti  =  {3} , 


(6.38a) 


(6.38b) 

(6.38c) 

(6.38d) 


(6.39a) 

(6.39b) 

(6.39c) 

(6.39d) 


while  the  right  structured  sub-matrices  can  be  defined  in  terms  of  Eq. 
then  the  corresponding  structured  sub-matrices  are: 


6.26 


If  we  write  V4  = 


Yi  y2  v3  y4 


y;  = 

V  1 

v2  v3  v4  ] 

,  (6.40a) 

r3  = 

v2 

v3  V4  ]  , 

(6.40b) 
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^  (N 
>11 

1 

V2 

(y2  +  v3)  v4 

,  (6.40c) 

Yi  = 

1 

V3 

(v2  +  v3  +  v4) 

(6.40d) 

Procedures  defined  above  create  matrices  that  are  orthogonal  in  both  dimensions,  but  normalizable  only  in  one. 
To  quantify  noise  resilience,  we  need  to  know  the  product  S  Calculating  is  trivial  because 

SS|4jj  is  diagonal,  while  calculating  is  more  challenging.  However,  for  4 ij  systems,  we  can  show  that 
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(6.42a) 

(6.42b) 

(6.42c) 

(6.42d) 


is  a  correct  inverse,  where 
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(6.43c) 

Noise  Resilience 

Following  the  same  additive  noise  model  as  prescribed  by  Eq.  |6.4|  we  can  project  the  noise  into  the  respective 
directions  of  the  sensor  spaced  V: 

n  =  W  +  ?n,  (6.44) 

where  the  pseudoinverse  can  be  written  within  the  context  of  4 ij  systems  as 


W  + 

=5^4  ij 


=  V 


:£  + 


Ns,4ij  ^a,4ij'M^s,4ij  ~  Xs>4u  =  Ys 


L  vs  4{j  contains  the  mapping  weights  of  information  for  each  of  the  vectors  of  V^: 


(6.45) 
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(6.46) 


Since  the  pMMP’s  sensor  space  contains  N  vectors,  L  v3  4i  also  contains  N  channels.  The  Euclidean  length  of 
each  of  those  vectors  represents  the  noise  magnitude  in  each  of  the  vectors  of  VSj4ij, 


Pv  =  \\l^ 


(6.47) 


For  each  measurement  set,  the  matrix  multiplications  reveal  that  each  of  the  vectors  making  up  Vs  4^  will  have 
easily  identifiable  noise  magnitudes: 


Py4  =  [  1  1  1  1]T,  (6.48a) 

P  v3  =  [  \/3  y/3  \/3  ] T  ,  (6.48b) 

P  Y2  =  [  1  V3]T,  (6.48c) 

r  j—  it 

Pvi  =  yf  •  (6.48d) 


Finally,  the  total  noise  magnitude  vector  is  the  concatenation  of  the  ones  defined  above, 


However,  since  the  intent  of  this  exercise  is  to  build  systems  that  perform  the  best  for  a  given  task,  it  also  follows 
that  it  would  be  desired  to  evaluate  system  performance  not  for  the  entire  sensor  space,  but  for  the  scene  space 
instead.'4  We  will  denote  that  space  y  to  match  its  computational  representation,  Y.  We  can  then  define  a 
transformation, 

X  =  X\Ys,4  ij,  (6.50) 

which  can  be  used  to  combine  N  measurements  into  the  vectors  approaching  y  or  estimating  Y, 


X  =  2¥W 


(6.51) 


Note  that  while  P  vs  4i  represents  the  noise  magnitude  for  the  reconstructable  vectors  represented  by  V^,  it 
would  be  incorrect  to  use  these  absolute  magnitudes  to  map  noise  from  reconstructables  to  the  desired  channels. 
Instead,  the  noise  characteristics  contained  within  Lv,  4i.  need  to  be  similarly  mapped: 
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(6.52) 


where  K  is  the  total  number  of  vectors  in  Y.  The  resulting  noise  magnitudes  within  those  vectors  can  be 
evaluated  in  a  philosophically  equivalent  way, 


Py  —  Illy  II2, 

and  then  can  be  combined  into  a  total  magnitude  vector 

T 

P.Y  =  [  Pyi  Py2  Pyjc  ]  • 


(6.53) 


(6.54) 


Space  Coverage 

To  properly  evaluate  a  given  partial  system,  it  is  important  to  know  not  only  the  system’s  noise  resilience,  but 
also  the  closeness  of  the  sensor  space  to  the  scene  space,  which  can  be  described  by  K  ordered  canonical  angles 
Ci  <  C2  <  •  •  •  <  Ck  1 '5  The  first  canonical  angle  Ci  is 

Cl  =  cos-1  (^  ^min^^  (xi  yl))  (6.55) 

Subsequent  canonical  angles  are  computed  by  evaluating  Eq.  |6.55|  with  the  portions  of  subspace  V  remaining 
after  the  elimination  of  vy.  The  best  case  scenario  is  when  <px  =  0,  which  means  JcV,  and  the  pMMP  spans 
the  desired  channels. 

While  Eq.  |6.55|  provides  an  intuitive  interpretation  of  the  canonical  angles,  there  are  more  efficient  ways  of 
computing  the  angles.  We  form  the  auxiliary  matrix 

1  =  X-Y(YTX)=  [  XI  x 2  •••  Xjf  ].  (6.56) 

and  compute  the  canonical  angles  from  the  singular  values  of  this  auxiliary  matrix  as 

Cfc  =  arcsin(crxj.  (6.57) 


6.4  Example  of  pMMP  Optimization 

To  find  the  best  pMMP  design  for  a  given  task,  we  need  to  optimize  for  both  noise  resilience  and  space  coverage. 
Because  those  properties  are  not  inherently  guaranteed  to  have  overlapping  minimums,  we  are  invariably  bound 
for  the  solution  to  be  a  point  on  the  Pareto  surface  of  a  multi-objective  optimization  problem.  We  have  found 
that  the  following  metric 

arg  min 

£ 


successfully  finds  appropriate  pMMP  designs.  The  choice  of  w,  {a^}  and  {/?*,}  provides  handles  to  adjust  the 
importance  of  all  the  various  parameters,  while  the  optimization  variable  vector  £  contains  six  values  to  construct 
three  generating  and  three  analyzing  vector  pairs.  The  first  four  variables  define  4>g11  6 (/>a1  and  9 a4  to 
produce  vectors  Gi  4-  and  Ai±,  while  the  second  two  variables  define  ipc  and  ip  a  to  prescribe  where  G2,±, 
A2  ±,  G3 and  A3 1±  reside  on  the  orthogonal  circles  with  respect  to  Gi  ±  and  Ai.±. 

To  illustrate  the  design  of  pMMPs,  we  consider  the  example  presented  in  Hoover  and  TyoJ19  Four  different 
coupons  of  an  ABS  plastic  material  were  exposed  to  different  fluences  of  high  energy  laser  energy,  and  the 
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(6.58) 


resulting  damaged  samples  had  their  monostatic  Mueller  matrices  measured  at  a  range  of  angles  from  —20°  to 
20°.  Performing  SVD  of  the  data  reveals  that  the  most  fundamental  three  measurement  channels  are 
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Note  that  the  original  paper  used  covariance  matrix  principal  component  analysis  that  resulted  in  a  different  set 
of  channels,  which  did  not  include  m 0o  in  any  of  the  measurements.  If  we  add  mOO  back,  then  the  maximum 
canonical  angle  between  the  two  spaces  is  3.1011°.  The  difference  is  small  enough  to  be  accounted  for  by  the 
extra  idealization  step  taken  in  Hoover  and  TyoP9 

We  used  MATLAB’s  built-in  genetic  algorithm  routine  together  with  Eq.  |6.58|  to  optimize  each  of  the  4 ij 
pMMP  designs  with  a*,  =  Pk  =  1  and  w  =  25.  Note  that  there  is  nothing  fundamental  about  our  choice  of  w 
—  we  tried  a  number  of  different  weights  between  1  and  100,  and  found  that  for  this  data  set  the  value  of  25 
provided  a  good  solution  where  the  space  coverage  penalty  was  just  significant  enough  for  the  reconstruction  of 
relevant  information  to  be  prioritized  over  the  noise  resilience. 

Table[4]shows  the  system  performances  that  we  were  able  to  find  for  each  of  the  defined  classes  of  polarimeters. 
We  can  point  out  that  the  space  coverage  seems  to  be  marginally  better  for  the  422  system  than  it  is  for  the  432 
or  the  433,  despite  the  latter  two  making  more  measurements  and  having  a  capacity  only  to  expand  the  space 
coverage  if  the  422  design  is  used  as  the  base.  That,  however,  is  purely  an  artifact  of  our  choice  of  w,  which 
leads  to  the  optimizer  finding  a  solution  with  slightly  better  noise  resilience  by  sacrificing  some  space  coverage. 
Practically,  the  designs  should  be  evaluated  on  whether  or  not  they  can  separate  the  different  objects  classes.  To 
determine  which  of  these  pMMPs  accomplish  that,  it  is  necessary  to  look  at  the  object  projections  onto  y.  We 
can  capture  this  by  looking  at  how  the  proximity  of  each  of  the  25  objects  from  each  of  the  four  types  of  objects 
to  the  nearby  classes  changes.  Instead  of  comparing  data  points  directly,  we  will  instead  piece-wise  interpolate 
the  comparison  class  and  determine  the  separation  for  each  object/class  both  in  y  and  y-. 
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where  a  and  /3  represent  the  object  classes,  7  represents  one  of  the  25  points  within  class  a,  and  8  represents 
one  of  the  24  line  segments  created  for  class  (3.  We  evaluate  the  geometric  mean  of  the  ratios  of  least  separation, 

1 

25 

(6.61) 


^o;,/3  — 


25 

n 

.7=1 


^  a, /3, 7, min 
^o:,/3,7,min 


When  ha,p  =  0,  classes  a  and  /3  have  collapsed  to  lie  on  top  of  each  other,  while  when  =  1,  the  separation 
between  classes  a  and  /?  has  remained  unchanged.  In  the  case  that  ha ^  >  1,  the  separation  within  the  recon¬ 
struction  is  greater  than  the  original  separation.  Although  this  presents  a  seemingly  interesting  scenario,  this 


Design 

N 

£ 1 

£2 

£l  -(-  W£2 

Qk 

hi, 2 

h.2,3 

^3,4 

hi, 3 

^2,4 

hi, 4 

400 

4 

1.911 

1.089 

29.142 

89.65° 

0.0003 

0.0006 

0.0006 

0.0001 

0.0001 

0.0000 

410 

5 

2.546 

0.717 

20.475 

56.90° 

0.1747 

0.2207 

0.2698 

0.3144 

0.1922 

0.0355 

411 

6 

2.669 

0.556 

16.569 

47.53° 

0.1770 

0.2417 

0.4435 

0.2232 

0.4506 

0.1438 

420 

6 

3.520 

0.469 

15.247 

43.18° 

0.5305 

0.4743 

0.5594 

0.2750 

0.5191 

0.5163 

421 

7 

3.528 

0.248 

9.727 

29.78° 

0.5184 

0.2960 

0.8056 

0.2087 

0.8787 

0.5441 

422 

8 

3.967 

0.002 

4.011 

2.32° 

1.0852 

0.9811 

0.9652 

0.9893 

0.9922 

0.9851 

430 

7 

3.318 

0.469 

15.042 

43.17° 

0.5133 

0.4793 

0.5625 

0.2705 

0.5435 

0.5151 

431 

8 

3.493 

0.249 

9.713 

29.84° 

0.5155 

0.3029 

0.8132 

0.1998 

0.8583 

0.5302 

432 

9 

3.932 

0.002 

3.987 

2.60° 

1.0940 

0.9921 

0.9632 

0.9893 

0.9914 

0.9869 

433 

10 

3.897 

0.002 

3.949 

2.57° 

1.1676 

0.9966 

0.9637 

0.9925 

1.0005 

0.9863 

Tabic  4:  Optimization  results  for  the  10  pMMP  system  classes.  The  optimization  targets,  e\  and  £2,  are  defined  in 
Eq.  |6.58]  while  Qk  represents  the  largest  canonical  angle.  The  values  for  ha ,p  are  calculated  via  Eq.  |6.6l1  Because  of  the 
way  that  the  four  classes  are  distributed  in  y,  knowledge  of  hi, 2,  h-2,3  and  hs,i  may  suffice. 


result  is  attributable  to  non-linearities  introduced  by  the  averaging  of  different,  space  projections  and  would  be 
compensated  by  another  ha ,p  elsewhere. 

Examining  the  performance  of  each  of  the  ten  pMMPs  optimized  designs  in  Table  [4] and  Fig.  |6.2|  it  becomes 
clear  that  the  422  system  is  the  first  design  of  the  defined  range  of  systems  that  accomplishes  the  task  of  matching 
the  space  coverage  and  thereby  separating  the  object  projections  adequately  for  object  detection. 

There  is  room  to  make  the  optimization  routine  more  elaborate.  For  example,  instead  of  matching  the  scene 
and  sensor  spaces  of  any  given  pMMP  class,  object  identification  can  be  done  in  the  measurement  space  itself. 
This  would  require  constructing  a  manifold  as  a  model  for  the  object  distribution  in  the  IV-dimensional  space, 
applying  the  proper  noise  model  and  looking  at  the  separability  of  the  classes  within  the  measurement  space. 
Performing  all  of  this  in  each  of  the  optimization  instantiations  is  computationally  intensive,  as  well  as  outside 
the  scope  of  this  development.  A  separate  discussion  is  warranted  to  address  that  level  of  optimization  properly. 

6.5  Conclusions 

Mueller  matrix  polarimeters  have  demonstrated  utility  recently  to  assist  in  target  identification,  and  the  use  of 
partial  Mueller  matrix  polarimeters  provides  a  way  to  develop  a  sensor  that  measures  the  polarization  featured 
needed  for  a  particular  detection  or  classification  task  without  having  to  measure  the  full  Mueller  matrix.  Previous 
designs  of  pMMPs  have  been  ad  hoc,  in  that  the  polarimeters  were  developed  by  hand.  In  some  instances,  there 
was  no  real  attention  paid  to  whether  or  not  the  pMMP  was  even  physically  realizable. 

We  have  developed  a  theory  of  pMMPs  that  enables  the  structure  of  a  pMMP  to  be  determined  from  the 
actual  generator/analyzer  pairs  used  to  form  its  instrument  matrix  W.  By  proper  analysis  of  W,  it  is  possible 
to  determine  the  portion  of  Mueller  matrix  space  that  a  particular  pMMP  measures.  We  developed  metrics  of 
optimality  for  pMMPs  that  are  based  on  balancing  their  SNR  performance  with  their  closeness  to  the  particular 
scene  space  at  hand.  The  performance  of  this  optimization  method  was  demonstrated  for  a  case  previously 
presented  in  the  literature.19 

7.  HYBRID  DOMAIN  MODULATED  POLARIMETERS 

Channeled  Mueller  matrix  polarimeters  and  the  concept  of  using  these  channels  was  first  introduced  by  Azzam. 16 
Azzam  published  a  very  specific  case,  1)  a  specific  temporal  framework  was  analyzed,  2)  an  implicit  assumption 
about  the  object  was  made,  the  object  had  no  temporal  bandwidth ,  i.e.,  the  object  was  stationary  in  time.  Oka, 
Sabatke,  Derniak,  Kudenov,  and  Hagen  then  demonstrated  both  spectrally  channeled  and  spatially  (over  spec¬ 
trum)  channeled  systems J!341]IH]lMll!2l  mostly  Stokes  polarimeters.  Dubreuil  et  al81  then  presented  a  spectrally 
channeled  Mueller  matrix  polarimeter,  which  of  course  is  non-imaging  since  the  focal  plane  array  is  used  to  resolve 
the  spectrum.  LaCasse,  Chipman,  Tyo,  and  LeMaster  and  HirakawcEEK]  then  described  micropolarizer  array 
partial  Stokes  polarimeters  as  channeled  systems,  and  LaCasse  et  al  presented  a  spatio-temporally  modulated 
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Figure  6.2:  Space  coverage  of  various  optimized  pMMP  designs.  Circles  represent  perfect  Mueller  matrix  object  projection, 
while  lines  are  the  approximation  that  each  pMMP  achieves.  Note  that  the  400  pMMP  collapses  all  measurements  onto  a 
single  line,  while  the  space  coverage  for  432  and  433  are  virtually  identical  to  that  of  422,  coupled  with  slightly  improved 
noise  resilience. 


hybrid  channeled  Stokes  system,'  and  subsequently  both  Myhre  et  and  Zhao  et  a/®  presented  spatially 
modulated  full  Stokes  polarimeters.  Finally  Alenin  and  Tyo1'  formalized  a  general  framework  which  describes 
channeled  polarimeters  almost  completely,  both  Mueller  and  Stokes. 

Prior  to  the  work  by  LaCasse  et  a/JZlE  10  bandwidth  in  channeled  polarimetric  systems  had  not  been  ad¬ 
dressed,  or  only  addressed  as  a  consequence  of  instrumental  “error."  Additionally,  prior  to  Alenin  and  Tyo1' 
channeled  systems  were  designed  in  an  ad-hoc  manner.  In  this  communication  we  address  bandwidth  using  the 
systematic  design  tools  introduced  by  Alenin  and  Tyo17  for  a  hybrid  spatio-temporally  modulated  channeled 
active  polarimetric  system. 

7.1  Formalism  and  channels 

Note  that  this  section  originally  appeared  in  our  other  publication  in  this  conference  proceeding 2Ji  and  is  de¬ 
rived/adapted  from  that  section  to  address  the  topic  of  this  communication  for  ease  of  reference.  Portions  may 
be  reproduced  verbatim,  however  quotes  will  not  be  used. 

We  use  the  Mueller-Stokes  mathematical  formalism  here,  as  it  is  most  commonly  used  in  instrumental  polar¬ 
ization  and  polarimeter  design.  This  analysis  is,  however,  agnostic  to  the  formalism  used,  a  coherence  formalisirPH 
with  periodic  modulators  could  also  be  used  and  would  have  similar  results.  In  the  next  sections,  it  should  be 
kept  in  mind  that  modulations  are  done  in  some  physical  domain,  they  are  periodic,  i.e.,  a  superposition  of 
sinusoidal  functions,  and  the  "channels"  are  the  resultant  ^-functions  which  ensue  from  the  Fourier  transform  of 
the  sinusoidal  modulations. 

7.1.1  Modulated  Mueller  formalism 

The  Stokes  parameters  are  described  by 


,  where  s0  >0,  Sq  >  s\  +  s\  +  Sg  (7.1) 


where  (•)  denotes  the  time  average,  so  is  proportional  to  the  total  irradiance,  si  is  proportional  to  the  preva¬ 
lence  of  horizontal  (0°)  over  vertical  (90°)  polarization,  s2  is  proportional  to  the  prevalence  of  +45°  over  —45° 
polarization,  and  s3  is  proportional  to  the  prevalence  of  right  circular  over  left  circular  polarization.^  Because 
optical  sensors  measure  a  quantity  proportional  to  the  time  averaged  Poynting  vector,  the  phase  information  is 
lost,  and  only  the  incoherent  time  averaged  polarization  information  can  be  obtained. 8 5 EH 

For  materials  which  can  be  described  via  linear  optical  interactions,  we  can  use  the  Mueller-Stokes  formalism. 
A  Mueller  matrix,  M,  is  a  matrix  which  linear  transforms  one  set  of  Stokes  parameters,  S;n,  into  another  set  of 
Stokes  parameters,  sout: 
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(7.2) 


Notice  that  M  £  M4x4  but  not  every  4x4  real  valued  matrix  is  a  Mueller  matrix  due  to  the  constraints  in 
Eqn]7.l|  see  GipH 

for  details. 


With  an  active,  or  Mueller  matrix,  polarimeteric  instrument,  we  must  modulate  in  irradiance  to  infer  the 
Mueller  matrix  of  an  object,  M(x),  where  x  =  [x  y  z  t  cr]  We  can  then  rewrite  Eqn. 


7.2 


to  have 

Mueller  matrices  and  Stokes  parameters  be  functions  of  space,  time,  and  wavelength  or  wavenumber.  '  Eqn|7.2| 
then  becomes 
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where  for  simplicity  we  fix  Sin.  Our  detector  then  measures  a  quantity  proportional  to  so,0ut(x).  For  a  Mueller 
matrix  measuring  instrument,  we  have  an  unknown  object  Mueller  matrix,  M  (x).  and  we  write  down  the 
instrument  equation  which  modulates  Stokes  parameters 


Sout(x)  =  A(x)  ■  MobjW  '  £(x) '  Sin 


(7.4) 


=  A(x) -Mobj(x)  •  S^.(x) 


(7.5) 


where  G(x),  A(x)  are  the  generator  and  analyzer  Mueller  matrices  respectively,  known  and  modulated  via  the 
physical  instrument.  The  generator  modulation  can  then  be  thought  of  as  only  a  Stokes  parameter  modulation, 
Sg(x). 


7.1.2  Channels 

Eqn  7.5  can  be  expanded  to  obtain  a  linear  equatioriPS for  So,0ut(x). 


3  3 

So, out  (x)  =^y^jaoi(x)sj(x)mij(x) 
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(7.6) 


where  aot(x)  are  the  elements  of  the  first  row  of  A(x),  Sj(x)  are  elements  of  sq(x),  and  rrqj(x)  are  elements  of 
Mobj(x).  We  can  then  take  the  Fourier  transform  of  So!0ut(x)  to  obtain 


3  3 

So, out  (p)  =  Aodp)  *  Sj(p)  *  Mij  (p) 
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(7.7) 


Where  x  — >  p  in  the  Fourier  transform,  *  denotes  convolution,  and  the  shift  to  capital  letters  indicates  a  function 
has  been  Fourier  transformed.  If  aot(x)  and  Sj(x)  are  superpositions  of  sinusoidal  functions,  then  Aoi(p)  *  Sj(p) 
is  a  set  of  5-functions,  and  each  Alij(p)  is  then  convolved  with  each  5-function  in  the  set.  The  complete  set  of 
5-functions  for  the  system 


3  3 

J2^2A0i(p)  *  Sj(p) 

2=0  j  —  0 

are  defined  as  the  channels  of  the  system,  or  the  system’s  channel  structured- 


(7.8) 


7.1.3  System  equation 

All  examples  in  this  communication  will  assume  a  quad-retarder  +  micropolarizer  array  Mueller  matrix  polarime- 
ter  system.  This  results  in  a  spatio-temporally  modulated  system  channel  structure.  The  details  of  the  design 
and  an  actual  instrument  implementation  are  in  Vaughn  et  al 24  and  the  system  equation  is  reproduced  below: 


s0ut  =E{x,y)  •  E(iT,  e4, 54)  •  g(v3,  e3, 53)  -M ohi(x,y,t)  •  S(^2, £2, S2)  •  E(iq,  ei,  5i) 

)sin  (7-9) 

where 

P(a:,  y )  =  micropolarizer  array  Mueller  matrix 

(7.10) 

R(z/j,  ej,  Sj )  =  retarder  Mueller  matrix 

(7.11) 

.  radians 

Vj  =  retarder  trequency  m  2 7 r - 

(7.12) 

ej  =  retarder  start  position  in  27t  radians 

(7.13) 

Sj  =  retarder  retardance  in  radians 

(7.14) 

M  (x,  y ,  t)  =  Mueller  matrix  of  the  object 

(7.15) 

Where  of  course  only  the  final  irradiance  value,  proportional  to  the  so,out(x)  element  of  sout(x) 

is  sampled.  The 

channels  are  encoded  in  so,out(x). 


7.2  Channel  design 

Designing  polarimetric  instruments  using  a  channeled  framework  is  somewhat  new  in  the  field,  especially  for 
instruments  which  are  not  spectrally  modulated.  Spectral  instrument  designers  have,  however,  utilized  channel 
design  in  an  iterative  way,  tweaking  the  system  designs  and  then  observing  the  channels  which  result,  then  again 
tweaking  the  system  design.  Here,  we  specify  constraints  and  needs,  and  then  optimize  the  system  directly  in  the 
channel  space  for  some  cost  function  dependent  on  those  specifications.  This  design  paradigm  allows  for  faster 
and  (sometimes)  conceptually  simpler  system  configuration. 

We  will  only  address  spatio-temporally  modulated  systems  here,  and  we  assume  modulations  are  of  the  type 

f(x,y,t)  =  h(x,y)  ■  g(t),  (7.16) 


i.e.,  that  any  modulation  is  mathematically  separable  between  time  and  space.  Non-separable  modulations 
can  be  constructed  (envision  a  rotating  focal  plane  array  with  a  micropolarizer,  or  a  spatial  light  modulator 
modulating  in  time  and  spatially),  but  we  will  address  these  in  a  future  publication.  A  separable  system  only 
allows  for  bandwidth  improvements  from  channel  cancellations  or  combinations,  while  a  non-separable  system 
may  allow  for  improvements  from  rotations  of  the  channel  structure,  however  the  latter  remains  an  open  question. 
We  emphasize  once  again  that  the  channels  are  (5-functions  in  the  Fourier  domain  which  result  from  sinusoidal 
modulations  of  irradiance  in  the  physical  domain ,  i.e.  space,  time,  wavelength,  etc. 


7.2.1  Notation 


Visualization  of  channel  structures  can  be  accomplished  by  graphing  the  sets  of  channels,  or  (5-functions  over  the 
Fourier  domain  dual  to  the  physical  modulation  domain.  Another  systematic  way  of  graphing  the  modulations 
is  the  frequency  phase  matrix  (FPM),  introduced  by  Alenin  and  Tyop  which  is  a  book  keeping  method  for 
the  channel  splitting  behavior,  the  signs,  and  real  and  imaginary  components  of  the  (5-functions.  In  this  com¬ 
munication  we  will  graph  the  channels  as  they  actually  appear  in  the  channel  space  to  reinforce  intuition  and 
understanding. 


positive 

negative 

0 

f— 1 

A 

V 

I  imag. 

> 

< 

Table  5:  Notation  for  channels. 


A  (5-function  can  be  characterized  by  its  position,  and  its  complex  magni¬ 
tude.  Table [hjoutlines  the  graphical  notation  that  we  will  use,  the  blue  triangles 
represent  the  real  part  of  the  magnitude,  red  triangles  represent  the  imaginary 
part  of  the  magnitude,  the  directions  that  the  triangles  point  represent  whether 
the  magnitude  is  positive  or  negative,  and  size  of  each  triangle  represents  the 
absolute  value  of  the  real  part  or  the  imaginary  part.  Note  that  we  will  only 
show  the  channels  for  a  single  Mueller  matrix  element  for  each  visualization, 
with  the  channels  for  all  other  Mueller  matrices  represented  by  light  gray-blue 
circles.  Keep  in  mind  that  many  channels,  for  each  Mueller  matrix  element, 
end  up  being  added  together  at  the  locations  shown.  Again,  the  relationship 
between  the  channels  of  different  Mueller  matrices  is  additive.  An  example  of 


the  channel  structure  for  a  spatio-temporally  modulated  system  for  77123  is  shown  in  Figj7.1| 


Some  of  the  examples  presented  here  will  be  normalized  to  a  temporal  frequency  range  of  [—1, 1],  this  is  because 
for  instrument  design  only  the  relative  frequencies  are  important,  we  are  always  limited  by  some  maximum 
sampling  rate  in  practice,  so  relative  bandwidth  with  respect  to  a  maximum  absolute  frequency  of  1  is  what 
must  be  optimized  for.  Due  to  the  assumption  of  separable  modulation  functions  for  space  and  time  and  Mueller 
matrix  physicality  conditions,  channels  are  fixed  to  travel  along  constrained  paths  in  the  Fourier  domain. 


7.2.2  Optimization 

Once  some  spatio-temporal  modulation  scheme  is  selected,  and  the  free  parameters  of  that  scheme  are  known, 
an  optimization  over  the  parameters  for  some  specific  cost  function  can  be  carried  out.  Alenin  and  TydlI1122J 
have  used  cost  functions  which  optimize  the  spectral  channel  structure  for  noise  performance,  but  other  cost 
functions  may  be  used.  In  this  communication,  we  jointly  optimize  for  bandwidth  and  noise  using  the  following 
cost  function: 


0{  P) 


[CN(p)]n 

dist(p) 


(7.17) 


£ 


Figure  7.1:  Example  of  a  spatio-temporal  channel  structure  with  5-functions  specific  to  77123.  The  maximum  bandwidth 
corresponds  to  the  minimum  distance  between  two  adjacent  channels,  taken  over  all  possible  adjacent  channel  pairs. 


where  p  is  a  vector  of  free  parameters  which  control  the  channel  structure,  the  dist(p)  function  quantifies  the 
distance  between  channels,  i.e.,  bandwidth,  CN( p)  is  the  condition  number  of  the  “mixing  matrix"  Q  (details 

about  Q  are  in  Alenin  and  Tyo1^),  and  n  £  R+  is  a  weighting  parameter  for  condition  number  corresponding  to 
noise  optimization.  This  is  valid  since  C N  >  1  by  definition.  The  optimization  then  minimizes  0( p)  over  p.  Note 
that  as  we  increase  n,  the  optimization  tends  to  favor  noise  performance  over  bandwidth.  A  system  conditioning 
metric  must  be  included  in  the  cost  function  to  ensure  reconstructablility  of  the  full  Mueller  matrix,  i.e.,  higher 
bandwidths  can  be  found  which  result  in  some  partial  Mueller  matrix  polarimeter  (pMMP)  reconstruction,  but 
the  system  conditioning  would  be  infinite  for  the  full  Mueller  matrix  reconstruction.  This  fact  has  utility  for 
pMMP  designs,  but  will  not  be  addressed  here. 


7.3  Bandwidth 

The  treatment  of  bandwidth  for  channeled  systems  is  mature  and  well  known  in  information  and  communications 
theory!®®1  The  difficulties  of  using  a  channeled  systems  framework  for  polarimetric  instruments  are  primarily 

1)  constructing  channels  in  2  or  more  dimensions  (many  systems  in  communications  theory  are  1  dimensional), 

2)  addressing  physicality  constraints  in  an  analytical  way,  and  3)  addressing  the  complicated  channel  mixing 
behavior  that  is  inherent  to  polarimetric  instruments.  3)  has  mostly  been  addressed  by  Alenin  and  Tyo.1'  2) 
is  a  complicated  subject  and  we  will  not  delve  into  details  here,  but  we  emphasize  again  that  these  constraints 
must  be  enforced  when  optimizing  for  some  cost  function. 

We  will  discuss  a  spatio-temporal  system,  with  modulation  in  the  domain 


(7.18) 


and  examples  will  be  from  a  quad-retarder  +  micropolarizer  array  system.24  First  we  will  briefly  review  convolu¬ 
tion.  Given  some  unknown  quantity,  m(f),  this  quantity  can  be  modulated  with  a  sinusoidal  function.  Without 
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Figure  7.2:  An  example  of  convolution  of  data  with  a  channel.  The  gray  band  represents  the  range  of  ec,  resulting  in  data 
being  outside  of  the  ec  range  for  —0.5  <  v  <  1.5.  This  implies  a  bandwidth  of  2  arb.  units  for  the  Mueller  data. 


loss  of  generality,  we  choose  cosine  here: 


f(t)  =  cos  27 rv0t  ■  m(t) 

F {/(*) =  F {cos2niy0t  ■  m(t)}t^u 
F{v)  =  F  {cos2tt *  F  {m(t)}t_>v 


F M  =  \ 


F{v)  = 


5(v  -  v0)  +  5(v  +  v0) 
M{v  —  v o)  M{v  +  v o) 


*  M(v) 


(7.19) 

(7.20) 

(7.21) 

(7.22) 

(7.23) 


where  .F{-}  is  the  Fourier  transform,  *  is  convolution,  and  6(v)  is  the  Dirac  delta  function.  The  last  line  is 
due  to  property  of  convolution  with  delta  functions.  This  gives  us  some  tools  for  conceptual  descriptions  for 
the  forward  problem  and  hence  the  bandwidth.  There  will  exist  a  set  of  channels  (^-functions)  for  each  Mueller 
matrix  element  Mij(p)  in  the  Fourier  domain.  For  each  channel  in  that  set,  My(p)  will  be  copied  at  that 
channel’s  location  with  MtJ  (0)  being  located  precisely  where  the  delta  function  is  located.  We  can  then  define 
the  bandwidth  of  for  some  threshold  ec  >  0  as  the  values  of  p  where  >  ec.  FigjD2]  clarifies  these 

concepts. 

The  polarimetric  system  channel  structure  contains  constraints  on  bandwidth  because  there  is  a  finite  dis¬ 
tance  between  channels  as  shown  in  FigjTT]  The  channel  structure  determines  the  bandwidth  available  for 
reconstruction.  When  the  bandwidth  of  the  data  becomes  greater  than  the  available  bandwidth,  then  channel 
crosstalk  occurs. 


7.3.1  Crosstalk 

Crosstalk  is  similar  to  aliasing,  but  not  the  same  phenomenon.  Crosstalk  is  the  result  of  limited  channel  band¬ 
width,  and  information  (convolutions  of  data)  in  the  channel  exceeding  the  bandwidth  of  that  channel  and 
"spilling  or  bleeding"  over  to  an  adjacent  channel.  Crosstalk  is  a  result  of  the  choice  of  channel  structure,  as 
opposed  to  the  sampling  rate  (aliasing),  even  a  continuously  sampled  channeled  (unaliased)  system  can  have 
crosstalk.  An  example  of  crosstalk  is  shown  in  Fig|7.3| 
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Figure  7.3:  An  example  of  channel  crosstalk.  Mueller  data  is  placed  at  two  channels,  with  the  distance  (bandwidth) 
between  them  less  than  the  bandwidth  of  the  Mueller  data.  When  added,  the  Mueller  data  from  different  channels 
adds  together,  leaving  no  remedy  to  differentiate  data  between  channels  in  the  region  of  bandwidth  crossover  when  given 
arbitrary  Mueller  data. 


7.3.2  Filtering 

Typically  crosstalk  can  be  “mitigated"  by  using  filters  around  the  channels  to  suppress  or  apodize  the  region 
where  crosstalk  occurs.  This  does  not  fully  mitigate  the  corruption  from  crosstalk  however  because 

•  Filters  which  apodize  in  some  way  result  in  smoothing  of  the  data,  essentially  removing  information. 

•  Similar  to  the  above,  apodization  or  cutoff  from  the  filters  essentially  reduces  the  bandwidth  of  the  resulting 
Mueller  data. 

•  Filters  won’t  help  much  in  the  case  where  a  great  deal  of  crosstalk  is  present. 

Filtering  is  needed,  but  cannot  fully  alleviate  the  crosstalk  issue.  Filtering  will  not  be  addressed  in  depth  here; 
the  literature  on  filtering  is  vast  and  mature  in  control  theory  and  electrical  engineering.  Keep  in  mind,  however, 
that  if  the  statistics  of  the  objects  being  measured  are  known,  then  optimal  filters  can  be  designed.^ 

7.3.3  Maximizing  bandwidth 

Our  focus  will  be  on  increasing  relative  bandwidth  to  reduce  crosstalk,  and  subsequently  increasing  the  system 
resolution  or  speed  for  spatio-temporally  modulated  active  polarimetric  systems.  In  order  to  maximize  the 
relative  bandwidth,  we  must  think  about  the  system  in  a  way  which  addresses  efficiency ,  otherwise  an  optimizer 
will  increase  the  maximum  frequency  (and  hence  the  relative  frequency  distance  between  channels)  ad  infinitum 
until  a  specification  is  met.  We  also  don’t  have  instruments  with  arbitrary  measurement  bandwidth.  In  order 
to  constrain  the  bandwidth  maximization  to  relative  frequencies,  we  can  normalize  all  of  the  channels  to  be 
contained  in  a  cube  (or  rectangular  prism  in  certain  cases)  where  the  maximum  frequency  is  normalized  to  be 
some  fixed  value.  We  can  choose  different  norms  to  accomplish  this  as  long  as  we  are  consistent.  The  two 
simplest  methods  are  1)  normalize  in  a  2-norm  way,  that  is  your  maximum  frequency  is  taken  as  a  vector  and 
normalized  by  its  2-norm  length,  and  all  other  channels  are  also  normalized  by  this  same  length,  or  2)  normalize 
in  an  oo-norm  way,  that  is  normalize  each  frequency  domain  coordinate  by  the  respective  maximum  frequency 
channel  domain  coordinate.  This  is  what  we  do  for  our  system  examples  here.  To  clarify  with  an  example, 


Figure  7.4:  An  example  of  varying  frequency  parameters  to  move  channels  around  in  the  channel  space.  Note  the 
cancellation/addition  at  certain  parameter  values  which  opens  up  more  bandwidth  between  the  channels.  Animated  in 
the  electronic  version,  use  the  controls  displayed  to  view. 


suppose  that  our  maximum  frequency  channel  is  located  at  [0.5,  0.5, 60]T,  then  the  normalization  factors  would 
be 

n2-„0rm  =  v7 0.52  +  0.52  +  602  (7.24) 


^oo-norm,0  =  0.5 
^oo-norm,  1  0.5 

^oo-norm,2  =  00 


and  for  an  arbitrary  channel  located  at  carb  =  [£arb  ?7arb  z'arb]  the  two  normalizations  would  be 

”  Car 
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(7.25) 

(7.26) 

(7.27) 
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In  the  examples  here,  we  only  normalize  the  locations  of  frequency  corresponding  to  the  temporal  domain,  v, 
because  the  examples  assume  a  fixed  micropolarizer  array  which  cannot  be  changed.24  This  results  in  normal¬ 
ization  of  the  spatial  frequency  coordinates  having  no  effect  on  the  analysis,  but  in  general,  if  optimization  over 
spatial  frequency  is  an  option,  the  spatial  frequency  channel  coordinates  would  also  need  to  be  normalized.  Nor¬ 
malization  ensures  that  for  a  relative  bandwidth  optimization  we  are  making  an  oranges  to  oranges  comparison 
as  channel  location  changes.  Note  that  the  use  of  different  normalization  types  will  lead  to  different  optimization 
outcomes. 


7.3.4  Channel  cancellation 

The  next  step  is  to  attempt  to  maximize  the  relative  channel  bandwidth,  now  that  the  channels  are  all  normalized 
to  a  maximum  frequency.  This  maximization  can  typically  be  accomplished  by  optimizing  over  the  system  channel 


structure’s  free  parameters  p  as  described  in  Eqn.  urn  Typical  parameters  for  a  spatio-temporally  modulated 
polarimeter  (which  is  not  spectrally  modulated)  include 

•  retardance  and  retarders,  spatially  or  temporally  modulated. 

•  spatial  or  temporal  analyzer/diattenuator  modulation. 

•  rotation  or  rotators;  these  elements  are  typically  combined  with  one  of  the  types  listed  above  to  achieve  a 
modulation. 

For  a  separable  channel  structure,  only  channel  cancellation/combination  or  reduction  of  overall  channels  may 
be  used  to  increase  the  relative  bandwidth.  Fig.  |7.4|  shows  channel  combination  as  relative  retarder  frequency 
is  changed  for  a  quad-retarder  +  micropolarizer  array  system.  At  certain  relative  frequencies,  channels  combine 
or  cancel  depending  on  their  magnitudes,  providing  larger  distance  (bandwidth)  between  channels.  Running  an 
optimizer  over  a  cost  function  can  then  find  an  optimal  channel  structure,  given  your  free  parameters  p.  An 
example  of  an  optimal  (to  the  best  of  our  knowledge)  channel  structure  for  the  quad  retarder  i§?  micropolarizer 
array  system  is  shown  in  Fig|7.5| 


.5 


Figure  7.5:  Optimal  channel  structure  for  a  specific  quad  retarder  +  micropolarizer  channel  structure  for  m23.  Reproduced 
from  a  figure  in  Vaughn  et  al.2'1 


7.3.5  Discussion 

Once  the  channeled  system  framework  is  understood,  and  the  free  parameters,  p  of  a  spatio-temporally  modulated 
system  are  specified  or  known,  then  it  is  straightforward  to  design  a  cost  function  and  run  an  optimizer  over  that 
function  to  optimize  for  bandwidth  or  jointly  for  bandwidth,  noise,  and  other  constraints.  The  most  difficult  part 
of  directly  optimizing  in  the  channel  space  is  not,  however,  defining  a  cost  function  and  running  an  optimization 
against  the  cost  function.  The  difficult  task  is  designing  a  model  which  properly  describes  the  channel  structure 
itself,  with  proper  physical  constraints.  We  have  designed  a  model  for  the  specific  case  of  a  quad-retarder  + 
micropolarizer  array  system,  but  we  hope  to  adapt  our  current  model  to  generate  generic  spatio-temporally 
modulated  systems  in  the  near  future. 

Additionally,  if  the  statistics  of  an  object  or  set  of  objects  are  known,  then  non-uniform  bandwidth  can  be 
maximized.  All  of  the  examples  shown  in  this  communication  optimize  for  an  equal  channel  bandwidth  between 
all  channels.  In  certain  cases,  more  bandwidth  may  be  wanted  for  certain  sets  of  channels  over  other  sets  of 
channels  and  for  certain  Mueller  matrix  elements.  This  can  all  be  accomplished  by  using  the  appropriate  cost 
function,  but  is  non-trivial  due  to  the  channel  mixing  which  occurs  between  Mueller  matrix  elements. 

7.4  Noise  and  sytematic  error 

The  channel  bandwidth  optimization  discussed  in  the  previous  section  must  be  paired  with  the  sensitivity  of  the 
system  to  deviations  from  the  ideal  optimized  parameters,  i.e.,  an  actual  system  will  have  non-ideal  modulations 
and  modulator  elements.  We  present  some  preliminary  results  on  systematic  deviations  here,  but  we  have  not 
fully  worked  out  how  to  measure  the  sensitivity  of  a  system’s  channel  structure  to  systematic  deviations  (often 
called  systematic  error  in  the  literature) .  The  sensitivity  to  random  noise  sources  is  also  discussed  here,  but  again 
we  have  not  worked  out  a  general  systematic  way  to  compute  the  sensitivity  of  channel  structures  to  random 
noise  when  bandwidth  is  taken  into  account  (including  filtering  effects).  Noise  effects  from  the  inversion  of  the 
Q  matrix,  however,  have  been  addressed.1' 

7.4.1  Systematic  error 

Systematic  error  will  occur  when  the  actual  system  deviates  from  the  designed  system,  or  when  there  is  some 
consistent  bias  due  to  physical  instrument  details.  Here  systematic  error  for  a  channeled  system  refers  to  the 
differences  between  the  real  system’s  channels  as  compared  with  the  ideal  channels  from  some  designed  channel 
structure.  Issues  arise  for  a  separable  channeled  system  when  channel  cancellation/combination  has  been  used 
as  a  tool  to  increase  system  bandwidth. 

For  separable  channeled  systems  a  potentially  serious  problem  arises:  at  the  locations  in  the  frequency  space 
where  the  channel  cancellation (s)  occurred,  channels  will  again  be  present  due  to  deviations  of  real  components. 
These  spurious  channels  will  be  convolved  with  Mueller  data,  and  introduce  channel  crosstalk.  We  really  only 
have  one  option  available,  to  minimize  the  magnitude  of  these  spurious  channels  so  that  the  crosstalk  into  the 
adjacent  channels  is  low.  An  example  of  is  shown  in  Fig.  |7.6|  the  small  triangles  represent  the  spurious  channels. 

We  have  found  a  way  to  reduce  the  effects  of  the  spurious  channels  for  the  quad  retarder  +  micropolarizer 
system,  which  was  to  re-optimize  over  our  channels  using  the  remaining  free  parameters  available  while  fixing 
the  ones  constrained  by  the  physical  instrument  components.  In  our  quad-retarder  +  micropolarizer  system 
example,  once  our  actual  retardances  were  fixed,  we  re-optimized  using  the  starting  position  of  each  retarder  and 
added  a  parameter  to  our  cost  function  which  characterized  the  magnitudes  of  the  spurious  channels  compared 
with  the  magnitudes  of  the  adjacent  channels.  For  general  spatio-temporal  systems  this  kind  of  method  will 
work  after  some  of  the  physical  components  are  specified,  if  there  are  any  free  parameters  left  to  optimize  over. 

In  fact,  an  iterative  approach  would  likely  yield  the  best  results  with  a  single  component  at  a  time  be¬ 
ing  sourced,  measured  and  characterized,  then  the  system  can  be  re-optimized  to  reset  the  other  components 
specifications  to  minimize  the  crosstalk,  and  then  the  process  is  repeated.8 
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Figure  7.6:  An  example  of  systematic  error  in  a  real  instrument,  the  small  triangles  represent  channels  which  are  present 
due  to  retardance  deviation  from  the  specifications. 
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Figure  7.7:  Gaussian  detector  noise  effects  on  channel  structure.  Note  that  these  are  the  same  channels  as  shown  in  Fig. 
|7.6|  but  plotted  along  the  v  lines.  The  left  column  is  the  real  part,  the  right  column  is  the  imaginary  part,  there  are  16 
colors  in  each  graph  representing  each  Mueller  matrix  element.  The  first  row  is  v  at  f  =  0,  r/  —  —0.5,  the  second  row  is  v 
at  £  =  —0.5,  r)  =  0  and  the  third  row  is  v  at  £  =  0,  r\  =  0.  Animated  in  the  electronic  version. 

7.4.2  Noise 

The  other  impacts  to  real  channeled  systems  are  random  noise,  for  instance  detector  noise.  The  results  presented 
here  are  preliminary,  and  we  only  modeled  detector  Johnson  (Gaussian  like)  noise  for  the  detector.  The  results 
are  simulated,  but  show  that  the  channel  structure  itself  is  quite  robust  to  detector  noise.  Here  we  switch  to  a 
different  view  of  the  channels  to  better  view  the  noise  effects.  Each  line  of  channels  in  the  v  direction  is  plotted 
in  Fig.  |7.7|and  Gaussian  noise  is  added  in  a  simulation  to  the  final  irradiance  for  this  result.  The  SNR  is  shown 
in  the  figure,  and  it  appears  that  the  channel  structure  is  stable  for  an  SNR  above  somewhere  between  1  and  2. 
This  result  will  need  to  be  validated  on  the  real  instrument,  and  we  plan  on  testing  this  in  the  near  future. 

We  have  not  modeled  the  effects  of  noise  on  the  reconstruction  for  our  specific  system  yet,  but  noise  effects 
have  been  addressed  for  channeled  systems  by  others EDEHEH 

7.5  Conclusion 

Channeled  polarimeter  design  has  changed  the  way  that  instrument  designers  approach  the  design  process,  and 
allows  engineers  to  systematically  design  both  general  and  task  specific  polarimeters.  We  have  presented  examples 
and  some  general  insight  into  using  channeled  system  design  for  spatio-temporally  modulated  Mueller  matrix 
polarimeters.  We  also  address  some  of  the  intricacies  of  channeled  design,  and  give  some  preliminary  examples 
of  possible  systematic  errors  and  noise  effects  on  channel  structures.  Furthermore,  we  give  an  example  of  how 
to  mitigate  some  of  the  systematic  errors,  and  in  the  future  will  use  a  different  approach  when  building  a  real 
instrument  by  re-optimizing  the  channel  structure  after  each  instrument  component  is  sourced.  In  the  future  we 
hope  to  build  a  fully  generic  channeled  system  model  for  spatio-temporally  modulated  instruments,  so  that  the 


community  can  design  their  own  instruments  in  a  systematic  way.  We  also  need  to  validate  our  systematic  error 
mitigation  on  our  real  instrument,  and  validate  the  effects  of  noise  on  the  channel  structure  of  a  real  instrument. 

8.  COHERENCE  MANIPULATION 

8.1  Polarization  and  its  Intrinsic  Relation  to  Coherence 

Current  design  philosophy  treats  conventional  imaging,  polarimetry,  spectrometry,  depth  imaging,  etc.,  as  sepa¬ 
rate  sensing  architectures.  However,  they  are  all  fundamentally  linked  through  the  statistical  description  of  the 
optical  field  embodied  in  the  mutual  coherence  function  at  the  sensor  aperture.  Once  the  optical  field  passes  the 
system  collection  aperture,  the  information  available  to  the  observer  has  been  fully  determined;  were  it  possible 
to  know  the  full  spatiotemporal  distribution  of  the  field  on  the  aperture,  then  the  output  of  any  conventional 
sensor  could  be  computed.  Unfortunately,  for  a  wide  class  of  sensing  systems  that  rely  on  incoherent  or  partially 
coherent  radiation,  this  complete  knowledge  is  impossible  to  obtain  at  optical  wavelengths  for  any  meaningful 
aperture  size.92  Since  current  optical  detectors  measure  second  moments  of  the  optical  field,  the  best  the  designer 
can  hope  to  accomplish  is  to  measure  some  meaningful  projection  of  the  field’s  mutual  coherence  function  across 
the  aperture. 

High  resolution  imagers,  spectrometers,  light  field  imagers,  polarimeters,  and  other  current  systems  use 
field  transducers  (lenses,  gratings,  polarizers,  etc.)  to  choose  the  projection  of  mutual  coherence  to  measure. 
Lens-based  imaging  systems  convert  spatial  Fourier  modes  on  the  aperture  into  irradiance  measurements  in  the 
focal  plane;  spectrometer  systems  convert  temporal  Fourier  modes  into  irradiance  measurements  on  the  detector 
plane;  extended  depth-of-field  imagers  use  aperture  coding  to  control  the  projections  of  the  light  field;  and 
various  field  modulation  strategies  allow  spectral,  polarization,  and  depth  information  to  be  demodulated  from 
channeled  intensity  data.  Recently  several  advanced  computational  and  compressive  sensing  strategies  that  are 
more  closely  tied  to  the  coherence  properties  have  emerged. 

A  key  aspect  of  our  research  program  is  study  of  the  theoretical  fabric  that  unifies  these  apparently  diverse 
sensing  strategies  through  the  theory  of  mutual  coherence.  Little  is  known  about  how  to  efficiently  measure 
specific  projections  of  the  field  coherence,  and  even  less  is  known  about  how  to  control  the  mutual  coherence 
properties  of  illumination  systems.  While  classical  statistical  optics  theory  describes  how  the  coherence  properties 
propagate  through  space,  the  tools  do  not  currently  exist  to  describe  the  transformation  of  coherence  upon 
interactions  after  scattering. 

Variable  coherence  tomography  (VCT)  was  developed  as  a  method  for  directly  measuring  the  second-order 
statistical  properties  of  quasi-homogeneous  mediaH?!94  VCT  uses  an  incident  wave  field  with  a  specific  mutual 
coherence  structure  that  allows  the  average  coherent  two-point  scattering  from  the  medium  to  be  probed.  VCT 
develops  a  speckle  field  that  is  uncorrelated  except  at  specific  lateral  offsets.  The  spatial  correlation  structure  of 
the  incident  beam  can  be  shaped  by  changing  the  parameters  of  the  source,  allowing  a  tomographic  reconstruction 
of  the  spatial  correlation  of  the  scattering  potential.  As  originally  developed,  VCT  is  based  on  scalar  scattering 
by  weak,  quasi-homogeneous,  isotropic  scattering  media.  However,  the  underlying  theory  is  equally  valid  for 
isotropic  and  anisotropic  scatterersP5  Variable  coherence  polarimetry  (VCPol)  was  introduced  by  our  group 
as  an  extension  of  the  theory  behind  VCT  to  include  polarimetric  information.16  VCPol  makes  it  possible  to 
probe  the  second-order  statistics  of  the  scattering  object  from  monostatic  measurements,  which  provides  greater 
capabilities  for  the  identification  and  classification  of  scene  abnormalities  in  remote  sensing  applications. 

8.1.1  Mutual  Coherence  Matrix 

The  second-order  correlation  properties  of  a  stochastic,  statistically  stationary  electromagnetic  beam  are  rep¬ 
resented  by  the  2x2  mutual  coherence  matrix.  This  matrix  is  the  electromagnetic  equivalent  to  the  mutual 
coherence  function  of  the  scalar  theory,^6  which  is  important  in  the  study  of  imaging  systems  with  partially  coher¬ 
ent  illumination.  The  generalized  van  Cittert-Zernike  theorem  is  a  fundamental  result  in  statistical  optics  that 
relates  the  spatial  coherence  properties  and  irradiance  distribution  of  a  beam  with  the  irradiance  distribution  and 
spatial  coherence  properties  of  the  quasi-homogeneous  source,  respectively,  via  a  Fourier  transform.  Equivalent 
results  in  the  context  of  the  joint  theory  of  coherence  and  polarization  have  been  proposed  recentlyP^SSl  The 
generalized  van  Cittert-Zernike  theorem  for  the  cross-spectral  density  matrix  developed  in  our  group  states  that 


Figure  8.1:  Diagram  of  a  source  designed  using  the  generalized  van  Cittert-Zernike  theorem  for  the  cross-spectral  density 
matrix  and  the  resulting  propagated  field. 


the  spatial  coherence  and  polarization  properties  in  the  far-field  of  a  quasi-homogeneous  electromagnetic  source 
are  given  by  the  Fourier  transform  of  the  functions  that  describe  its  polarization  and  spatial  coherence  properties. 
The  importance  of  this  theorem  in  both  the  scalar  and  vectorial  theories  is  that  it  allows  the  design  of  sources 
with  prescribed  spatial  coherence  distributions.  In  the  case  of  the  vectorial  theory,  this  capability  is  further 
extended  to  the  polarization  state  of  the  beam.^  Therefore,  knowledge  of  the  propagation  relations  followed 
by  the  mutual  coherence  matrix  is  fundamental  in  the  design  of  sources  with  prescribed  spatial  coherence  and 
polarization  properties  that  can  be  used  in  the  implementation  of  the  VCPol  technique.  Fig.  |8.1|  provides  an 
example  of  a  source  with  such  characteristics.  The  cross-spectral  density  matrix  depicted  in  the  figure  represents 
a  beam  that  is  unpolarized  in  the  usual  one-point  sense  but  polarized  in  the  two-point  sensed 

8.1.2  Statistical  description  of  scattering 

The  study  of  the  interaction  between  an  incident  field  and  a  random  object  is  often  carried  out  in  a  deterministic 
way  for  a  pair  of  single  realizations  of  the  field  and  object  and  then  repeated  for  a  large  number  of  different 
realizations.  In  this  Monte  Carlo  approach,  the  final  result  is  obtained  as  the  ensemble  average  of  the  results 
obtained  for  the  different  realizations.  An  alternative  strategy  is  to  describe  both  the  field  and  the  random  object 
as  ensemble  averages  from  the  beginning  and  find  the  scattering  relation  that  they  satisfy.  A  first  approximation 
to  the  description  of  this  interaction  between  a  stochastic,  statistically  stationary  field  and  a  random  object,  based 
on  the  first  Born  approximation,  can  be  found  in  the  literature  for  the  scalar  caseJHtt23ll££l  The  corresponding 
result  for  the  electromagnetic  case  has  been  addressed  in  earlier  work  by  our  group.16  This  approach  to  the 
calculation  of  the  scattered  field  in  statistical  optics  is  more  computationally  convenient  than  the  Monte  Carlo 
method  because  it  uses  the  statistical  properties  of  the  random  object  and  stochastic  incident  field  to  predict 
the  scattered  field  directly.  The  main  inconvenience  of  this  approach  is  that  the  scattering  potential  and  cross- 
spectral  density  matrix  of  the  incident  field  must  have  a  relatively  simple  mathematical  representation  in  order 
to  compute  the  scattered  field  in  closed  form.  Examples  of  such  relatively  simple  representations  are  the  quasi- 
homogeneous  fields  and  media,  which  constitute  a  good  approximation  in  many  applications!!®01  To  account 
for  situations  in  which  the  first  Born  approximation  is  not  applicable  (e.g.,  multiple  scattering  and  scattering 
in  the  near-field)  the  scattering  theory  must  be  extended.  The  extension  of  the  statistical  scattering  theory  to 
cover  these  situations  is  part  of  the  work  that  will  be  done  in  future  research. 

We  have  developed  a  design  for  a  novel  active  sensing  system  able  to  measure  the  elements  of  the  mutual 
coherence  matrix,  allowing  us  to  obtain  a  proof  of  concept  of  VCPol.  Fig.  |8.2|is  an  example  of  a  sensing  systems 
able  to  control  the  second  order  statistical  properties  of  the  illumination  and  measure  the  elements  of  the  mutual 
coherence  matrix  resulting  from  its  interaction  with  a  specimen.  The  system  combines  a  subsystem  to  control 
the  second  order  statistical  properties  of  the  illumination  with  a  Stokes  polarimeter  formed  by  variable  retarders 
VR1  and  VR2,  and  linear  polarizer  P4.  Measurements  with  different  second-order  correlation  properties  of  the 
source  may  be  used  to  compute  the  elements  of  the  mutual  coherence  matrix.  The  system  in  Fig.  |8.2|  is  not 
the  only  configuration  able  to  measure  these  quantities  and  we  must  examine  other  possibilities  to  determine  an 
accurate  and  efficient  way  to  measure  the  elements  of  the  mutual  coherence  matrix.  Furthermore,  an  analysis 
of  the  amount  of  data  required  and  the  development  of  the  methods  necessary  to  reconstruct  the  VCPol  results 
from  measurements  of  the  mutual  coherence  matrix  are  two  of  our  ongoing  and  future  research  objectives. 


Figure  8.2:  Diagram  of  a  novel  sensing  system  to  indirectly  measure  the  elements  of  the  mutual  coherence  matrix. 
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